A unifying mechanism governing interbrain neural relationship during social interactions
Abstract
A key goal of social neuroscience is to understand the interbrain neural relationship—the relationship between the neural activity of socially interacting individuals. Decades of research investigating this relationship have focused on the similarity in neural activity across brains. Here, we instead asked how neural activity differs between brains, and how that difference evolves alongside activity patterns shared between brains. Applying this framework to bats engaged in spontaneous social interactions revealed two complementary phenomena characterizing the interbrain neural relationship: fast fluctuations of activity difference across brains unfolding in parallel with slow activity covariation across brains. A model reproduced these observations and generated multiple predictions that we confirmed using experimental data involving pairs of bats and a larger social group of bats. The model suggests that a simple computational mechanism involving positive and negative feedback could explain diverse experimental observations regarding the interbrain neural relationship.
Introduction
What is the relationship between the neural activity of socially interacting individuals? This central question in social neuroscience has motivated nearly two decades of research spanning a diversity of species and methodologies (e.g. Babiloni and Astolfi, 2014; Dumas et al., 2011; Freiwald, 2020; Hasson et al., 2012; Hasson and Frith, 2016; Hoffmann et al., 2019; Kingsbury and Hong, 2020; Koike et al., 2015; Konvalinka and Roepstorff, 2012; Liu et al., 2018; Montague et al., 2002; Redcay and Schilbach, 2019; Scholkmann et al., 2013; Schoot et al., 2016; Testard et al., 2021; Tseng et al., 2018; Wass et al., 2020). Yet, despite this diversity, nearly all research has tackled the study of interbrain relationship from a single perspective: considering the neural activity of two interacting individuals as the two variables of interest and searching for similarities between them. Commonly, similarity was assessed using measures related to either correlation (Dikker et al., 2014; Kawasaki et al., 2013; KingCasas et al., 2005; Kingsbury et al., 2019; Kinreich et al., 2017; Levy et al., 2017; Liu et al., 2017; Montague et al., 2002; Piazza et al., 2020; Silbert et al., 2014; Spiegelhalder et al., 2014; Stephens et al., 2010; Stolk et al., 2014; Tomlin et al., 2006; Zadbood et al., 2017; Zhang and Yartsev, 2019) or coherence (Cui et al., 2012; Dikker et al., 2017; Dumas et al., 2010; Goldstein et al., 2018; Levy et al., 2017; Lindenberger et al., 2009; Montague et al., 2002; Mu et al., 2017; Stolk et al., 2014; Yang et al., 2020; Yun et al., 2012).
However, one aspect of the interbrain neural relationship has remained unexplored: the difference in neural activity between brains. We reasoned that interbrain difference is more than simply a lack of interbrain correlation. Specifically, the detailed dynamics of the difference between brains could contain information not captured by measures of interbrain similarity such as correlation. Therefore, here we took an approach that focused on both the similarity and difference in interbrain activity, and on how the two coevolve over time. We applied this approach to neural activity simultaneously recorded from socially interacting Egyptian fruit bats (Rousettus aegyptiacus), a mammalian species known for its high level of sociality (HerzigStraschil and Robinson, 1978; Harten et al., 2018; Prat et al., 2015; Prat et al., 2016; Prat et al., 2017; Omer et al., 2018; Cvikel et al., 2015; EgertBerg et al., 2018; Kwiecinski and Griffiths, 1999). Our recordings targeted the frontal cortex (Figure 1—figure supplement 1), a region implicated in social cognition in rodents, bats, nonhuman primates, and humans (Adolphs, 2001; Amodio and Frith, 2006; Cao et al., 2018; Chang et al., 2013; Eliades and Miller, 2017; Forbes and Grafman, 2010; Haroush and Williams, 2015; Kingsbury et al., 2019; Liang et al., 2018; Miller et al., 2015; Nummela et al., 2017; Ong et al., 2021; Pearson et al., 2014; Rose et al., 2021; Rudebeck et al., 2008; Tremblay et al., 2017; Zhang and Yartsev, 2019; Zhou et al., 2017). We then used modeling to better understand the experimentally observed interbrain relationship, which in turn produced a number of predictions that we subsequently confirmed using experimental data. Combined, the insights from the experimental and modeling approaches provide a unifying explanation to a range of experimental observations regarding the interbrain relationship during social interactions.
Results
Interbrain relationships in social and nonsocial contexts
We performed wireless extracellular neural recording simultaneously from pairs of bats. We recorded four types of neural signals: local field potential (LFP) power in the 30–150 Hz and 1–29 Hz bands (previously identified as relevant frequency bands in bat frontal cortical LFP [Zhang and Yartsev, 2019]), multiunit activity, and single unit activity (Materials and methods). The analyses in the main text focus on 30–150 Hz LFP power, as previous work has shown that this frequency range exhibits strong interbrain correlation in bats (Zhang and Yartsev, 2019).
To compare the interbrain relationship between social and nonsocial conditions, neural activity was recorded in two sets of experiments. In one experiment, pairs of bats behaved freely and interacted with each other inside a chamber (‘onechamber sessions’; Figure 1A). In the second set of experiments, the same bats freely behaved in separate, identical chambers (‘twochambers sessions’; Figure 1B). There were three types of twochambers sessions: (1) two bats each freely behaving in isolation; (2) two bats each freely behaving in the presence of identical auditory stimuli (playback of bat calls); (3) two bats each freely behaving and interacting with a different partner in separate chambers.
Plotting the neural activity of the two bats in onechamber sessions shows a high degree of similarity between brains (Figure 1C), which is not the case for twochambers sessions (Figure 1D), as demonstrated previously (Zhang and Yartsev, 2019). Such plots highlight the degree of interbrain similarity, but makes it easy to overlook the detailed dynamics of the difference between brains. We therefore sought an analysis framework that would enable us to explicitly examine interbrain difference and similarity side by side.
Relative magnitudes and timescales of interbrain difference and mean components
When studying the interbrain relationship, two typical variables of interest are the neural activity of each brain (e.g. in Figure 1C, each variable is the normalized LFP power averaged across the recording channels of one brain). We can represent them as a twodimensional vector $\left(\begin{array}{c}{\mathit{a}}_{1}\left(t\right)\\ {\mathit{a}}_{2}\left(t\right)\end{array}\right)$ , where ${\mathit{a}}_{1}\left(t\right)$ and ${\mathit{a}}_{2}\left(t\right)$ are the activity of bat 1 and bat 2 at time $t$, respectively. Through a change of basis, the same activity can be represented under another orthogonal basis as the mean and difference between the two brains: $\left(\begin{array}{c}{\mathit{a}}_{M}\left(t\right)\\ {\mathit{a}}_{D}\left(t\right)\end{array}\right)=\frac{1}{2}\left(\begin{array}{c}{\mathit{a}}_{1}\left(t\right)+{\mathit{a}}_{2}\left(t\right)\\ {\mathit{a}}_{1}\left(t\right){\mathit{a}}_{2}\left(t\right)\end{array}\right)$ , where ${\mathit{a}}_{M}\left(t\right)$ is the mean component of the activity, and ${\mathit{a}}_{D}\left(t\right)$ is the difference component (here the difference component is defined as $\frac{1}{2}\left[{\mathit{a}}_{1}\left(t\right){\mathit{a}}_{2}\left(t\right)\right]$ rather than ${\mathit{a}}_{1}\left(t\right){\mathit{a}}_{2}\left(t\right)$ so as to have the same scale as the mean component). ${\mathit{a}}_{M}\left(t\right)$ represents the common activity pattern shared between the brains, that is, what is similar between the brains; ${\mathit{a}}_{D}\left(t\right)$ represents the momenttomoment activity difference between the brains. Note that the two components can vary independently—a large (or small) ${\mathbf{a}}_{M}\left(t\right)$ does not necessarily imply a small (or large) ${\mathbf{a}}_{D}\left(t\right)$ .
Figure 1E shows the mean and difference activity components for the example session from Figure 1C. Visualizing the two components sidebyside, two relationships are immediately apparent. First, the mean component had much larger variance than the difference component. This is expected given the high degree of similarity in neural activity between the brains: while the variances are not entirely determined by interbrain correlation, a positive correlation does mathematically imply larger variance for the mean compared to the difference component (see Materials and methods for details). Second, the difference component evolved over time at much faster timescales than the mean component. This is even more apparent when examining the autocorrelations of the mean and difference components (Figure 1G), where the narrower autocorrelation of the difference indicates that it varied faster than the mean.
In contrast, the twochambers sessions showed a very different picture. Figure 1F shows the mean and difference components for the example session from Figure 1D. Here, the variances of the two components were comparable, and so were their timescales, which is also apparent from their comparable autocorrelations shown in Figure 1H.
We next quantified the magnitudes and timescales of the mean and difference components for 30–150 Hz LFP power on all sessions, using variance as a measure for magnitude, and power spectral centroid as a measure for timescale (Figure 1I–L). The power spectral centroid was calculated as follows. Given the time series of an activity component, we subtracted from it its average over time before calculating its power spectrum, and then computed the weighted average of the frequencies, where each frequency was weighted by the power at that frequency (Materials and methods). The relationships seen in the onechamber example above was robust: on every onechamber session, the difference component had smaller magnitudes (Figure 1I) and faster timescales (Figure 1J) than the mean component (note that the relative magnitudes are a reflection of interbrain correlation). On twochambers sessions, in contrast, the mean and difference were comparable (Figure 1K–L). These relationships are summarized in Figure 1M–T for all four neural signals: LFP power in the 30–150 Hz and 1–29 Hz bands, multiunits, and single units all showed the same significant trends (see Figure 1—figure supplements 2–4 for examples and detailed results for 1–29 Hz LFP power, multiunits, and single units).
Thus, the interbrain neural relationship during social interactions is characterized by two robust signatures of the mean and difference components: their relative magnitudes (which reflects interbrain correlation) and their relative timescales. As mentioned above, the observed relative magnitudes are mathematically implied by the presence of interbrain correlation (Materials and methods). It is important to determine whether the observed relative timescales are also implied by interbrain correlation and the observed relative magnitudes. As we show in Materials and methods (section ‘Surrogate data and the relationship between interbrain correlation and mean and difference components’), the answer is no: having a given combination of correlation, mean component variance, and difference component variance does not place constraints on the timescales of the two components. Specifically, it does not constrain the difference to be faster than the mean. To explicitly demonstrate this, we constructed a set of surrogate data as a counterexample (Figure 2). The surrogate data was designed to have interbrain correlation, mean component variance, and difference component variance that are identical to the actual data. However, unlike the actual data, the surrogate data was designed to have difference components that were slower than the mean components, the opposite of what was experimentally observed (Figure 2C and E). Thus, the relative timescales of the difference and mean components are not dictated by their relative magnitudes or by interbrain correlation.
What, then, might explain the robust relationships observed between the timescales and magnitudes? And are there separate mechanisms responsible for the observed relative timescales on the one hand, and the observed relative magnitudes (and the related phenomenon of interbrain correlation) on the other? To address these questions, we next model the observed neural activity to infer the computational mechanisms governing the interbrain difference and mean components, and by extension, mechanisms underlying interbrain correlation.
Model suggests a feedback mechanism governing interbrain neural relationship
We modeled the neural activity of two bats using a linear differential equation (see Figure 3A for the equation; Dayan and Abbott, 2005). In the model, the neural activity variables interact through the functional coupling matrix $\mathit{C}=\left(\begin{array}{cc}{\mathit{C}}_{S}& {\mathit{C}}_{I}\\ {\mathit{C}}_{I}& {\mathit{C}}_{S}\end{array}\right)$, where $\mathit{C}}_{S$ is the strength of functional selfcoupling and $\mathit{C}}_{I$ is the strength of functional acrossbrain coupling. When modeling onechamber sessions, ${\mathit{C}}_{I}>0$, so that the activity of each bat influences the other bat’s activity. This functional acrossbrain coupling models the effects of sensorimotor interactions and attentional processes (Hasson et al., 2012): for example, when bat 1’s neural activity increases due to its active movements (Gervasoni et al., 2004; McGinley et al., 2015), the movements create sensory inputs to bat 2, which can drive bat 2’s neural activity to the extent that bat 2 is paying attention (Chun et al., 2011; Driver, 2001; Fritz et al., 2007; Reynolds and Chelazzi, 2004). On the other hand, when modeling twochambers sessions, ${\mathit{C}}_{I}=0$, so that the activity of each bat does not influence the other bat’s activity. Furthermore, in both onechamber and twochambers models, the activity of each bat is modulated by its own behavior, which is simulated using Markov chains based on empirical behavioral transition frequencies (Materials and methods). Note that the usage of the term ‘functional acrossbrain coupling’ in our model should be distinguished from the sense it is sometimes used in the literature to simply denote a presence of neural correlation or coherence across brains (e.g. Levy et al., 2017).
This model was used to simulate neural activity, which was then analyzed using the same methods used to analyze experimental data. The model qualitatively reproduces the experimentally observed neural activity patterns: simulated activity shows variability over time, modulation by behavior, and correlation across brains (Figure 3B; see Figure 3—figure supplement 1 for a quantitative comparison between model and data). More importantly, the model also qualitatively reproduces all the experimentally observed relationships between the magnitudes and timescales of the difference and mean components, as can be seen in the activity and autocorrelation of the two components on an example simulation (Figure 3C–D). Specifically, the difference component had smaller magnitudes and faster timescales than the mean component in the onechamber model, but not in the twochambers model (Figure 3C–F).
What mechanisms in the model are responsible for these results? We now analyze the model to answer this (see Materials and methods for the detailed analysis). The model describes the evolving neural activity of two bats and uses separate variables to represent the activity of each bat. This is a basis in which the neural activity variables are coupled to each other in the onechamber model: the activity of bat 1 influences the activity of bat 2, which in turn feeds back on the activity of bat 1. It is much easier to understand how activity evolves if we change to a basis in which the activity variables are uncoupled. This basis consists of the eigenvectors of the functional coupling matrix $\mathit{C}$. Under the eigenvector basis, the neural activity variables become the mean activity across bats, and the difference in activity between bats (Figure 3G). Thus, the uncoupled activity variables are precisely our variables of interest: the interbrain mean and difference components. Each of the two components provides feedback onto itself, with feedback strengths being the eigenvalues of $\mathit{C}$: $\left({\mathit{C}}_{S}{\mathit{C}}_{I}\right)$ for the mean component and $\left({\mathit{C}}_{S}+{\mathit{C}}_{I}\right)$ for the difference component. In the twochambers model, ${\mathit{C}}_{I}=0$, so the two components receive equal feedback. In the onechamber model, ${\mathit{C}}_{I}>0$, so functional acrossbrain coupling acts as positive feedback for the mean component, which amplifies the mean component while slowing it down. On the other hand, functional acrossbrain coupling acts as negative feedback for the difference component, which suppresses the difference component while speeding it up (Figure 3H). Thus, in the model, a single mechanism—opposite feedback, that is, positive feedback to the mean component and negative feedback to the difference component—contributes to all of our observations relating the magnitudes and timescales of the mean and difference components.
We can gain a more precise understanding of the dependence of model output on functional coupling, by analyzing a reduced version of the model. The reduced model assumes a simplified structure for behavioral modulation, which allows derivation of simple analytic expressions for the variance ratio and power spectral centroid ratio of the mean and difference components (Materials and methods). Specifically, the variance ratio and centroid ratio can be shown to be approximately $\frac{{\mathit{C}}_{S}+{\mathit{C}}_{I}}{{\mathit{C}}_{S}{\mathit{C}}_{I}}$ and $\frac{{\mathit{C}}_{S}{\mathit{C}}_{I}}{{\mathit{C}}_{S}+{\mathit{C}}_{I}}$, respectively (Figure 3I–J). This simple dependence on the functional coupling parameters $\mathit{C}}_{S$ and $\mathit{C}}_{I$ is consistent with our analysis of the full model above: functional acrossbrain coupling $\mathit{C}}_{I$ modulates the mean and difference components in opposite directions. Furthermore, it identifies the parameter regime in which the experimental observations lie—the parameter regime of positive $\mathit{C}}_{I$ (Figure 3I–J).
In the full version of our model, in addition to functional coupling, behavioral modulation also contributes to the interbrain relationship (the effects of behavioral modulation are analyzed in depth in Materials and methods and Figure 3—figure supplement 2). In particular, the tendency of bats to engage in the same behavior at the same time (Figure 3—figure supplement 2G) contributes to interbrain correlation. This led us to ask whether this form of behavioral coordination alone is sufficient to explain experimental data in the absence of functional acrossbrain coupling. One specific experimental observation suggests that the answer is no: Figure 4A, B and G show that, after removing time periods of coordinated behavior from experimental sessions, interbrain correlation persisted. We found that our model can reproduce this phenomenon (Figure 4C, D and G). Next, we simulated an alternative model with behavioral coordination, but without functional acrossbrain coupling (Materials and methods). Figure 4E–G show that, in the absence of functional acrossbrain coupling, interbrain correlation disappeared after removing time periods of coordinated behaviors. Thus, in the framework of our model, functional acrossbrain coupling is necessary to reproduce the experimental results.
In our model, functional acrossbrain coupling took a simple form: the activity variables of the two brains are linearly coupled. We made this modeling choice for its simplicity, but one may wonder, does the form of functional acrossbrain coupling matter, or is any form of coupling capable of reproducing the experimental results? While we cannot exhaustively test all alternative ways of modeling functional coupling, we sought to explore one alternative, the Kuramoto model, because it is a wellknown model of coupling between oscillatory variables (Strogatz, 2000; Acebrón et al., 2005). The Kuramoto model has been widely used to model diverse synchronization phenomena in physics, chemistry, and biology, including neural synchronization both within and across brains (e.g. Breakspear et al., 2010; Cumin and Unsworth, 2007; Dumas et al., 2012; Schmidt et al., 2015). Here, we adapted it to model interbrain synchronization in bats. In this model, the fluctuating neural activity of the interacting bats are abstracted as oscillators whose phases are dynamically coupled depending on their phase difference (Materials and methods). This phasecoupling mechanism is able to reproduce interbrain correlation and the relative magnitudes of the mean and difference components from the data; importantly, however, it does not reproduce the relative timescales of the mean and difference components from the data (Figure 4—figure supplement 1). This demonstrates that not all forms of functional coupling are equally able of reproducing the interbrain relationship.
In summary, our main model provides a simple explanation for a set of robust, yet puzzling relationships between the interbrain mean and difference components. Namely, through opposite feedback to the mean and difference components, functional acrossbrain coupling simultaneously modulates both their magnitude and timescales in opposite directions—the amplification and slowing down of the mean and the suppression and speeding up of the difference are all manifestations of a single mechanism.
Testing predictions emerging from the model
In addition to explaining experimental observations, the model also makes novel predictions about previously unexamined aspects of the interbrain relationship. In this section, we describe three predictions about disparate aspects of the interbrain relationship: (1) correlation between neural activity variables defined by rotations of the twobat activity space axes; (2) relationship between the difference and mean components across sessions, rather than within sessions; (3) the encoding of the behaviors of one bat by the neural activity of the other bat. While the three predictions are seemingly unrelated to one another, they are in fact all motivated by the same model mechanism of functional acrossbrain coupling. After describing each prediction, we turn back to the data to test it. If the data can confirm the predictions, it would support the validity of the model.
In the model, the neural activity of each bat is a variable. These two activity variables can be thought of as evolving within a twodimensional space, that is, a space whose axes are the neural activity of each of the two bats (Figure 3G left). The evolution of activity in this space can be equivalently described using alternative sets of activity variables: for example, by rotating the axes by 45°, the activity variables become the difference and mean components (Figure 3G right). Rotating the axes in this way changes not only the activity variables, but also the strength of functional coupling between the variables. In the onechamber model, if we rotate the axes smoothly, the strength of functional coupling also changes smoothly (Figure 5A). In particular, the coupling strength is positive and at its maximum when the axes correspond to the activity of each bat, and it decreases to zero as the axes rotate by 45° (corresponding to the mean and difference components), and finally becoming negative as the axes rotate further (Figure 5A). After any amount of rotation, we can calculate the correlation between the activity variables. The model predicts that, as the axes rotate, the correlation between the neural activity variables would mirror the coupling strength (Figure 5B). This effect is not due to the changing behavioral modulation upon axes rotation, as can be seen after regressing out the influence of behavior from the activity (brown curve in Figure 5B; see Materials and methods for procedures used to regress out behavior). As we turn back to the data, we found that it clearly confirmed the model predictions: correlation as a function of rotation showed the same relationship as in the model (Figure 5C). Similarly, data from the twochambers sessions (Figure 5F) also confirmed the predictions of the twochambers model (Figure 5D–E): namely, near zero correlation across all axes rotation angles.
In the previous sections, we have examined in detail the relationships between the magnitudes and timescales of the interbrain activity components in single sessions, but the model also makes predictions regarding their relationship across sessions. As described above and in Materials and methods, functional acrossbrain coupling acts as positive feedback for the mean component and negative feedback for the difference component, modulating both their magnitudes and timescales in opposite directions (Figure 3H). Thus, the model predicts that the relative magnitudes and relative timescales of the two components are tied together by this mechanism and do not vary independently. This can be seen by plotting the relative magnitudes against the relative timescales for different simulations, which shows linear relationships between them (Figure 5G; see Materials and methods for detailed analysis of these relationships). Turning to the data, we see similar linear relationships, again confirming the predictions of the model (Figure 5H). It is important to note that these are very specific predictions that were not at all obvious or expected from our previous knowledge of the interbrain relationship. The fact that they were verified by the data provides strong support for the validity of the model.
In the model, the neural activity of each bat is directly modulated by its own behavior; moreover, it is indirectly modulated by the behavior of the other bat, through functional acrossbrain coupling. This suggests that the neural activity of each bat should represent the behavior of the other bat independently from encoding its own behavior. To quantify this in the model, we first regressed out the behavior of each bat from its own neural activity, then examined to what extent the residual neural activity encodes the behavior of the other bat, using the receiver operating characteristic (ROC) curve. The reason we first regress out each bat’s behavior is the following. If two bats happened to engage in the same behavior at the same times, the activity of each bat would appear to encode the behavior of the other bat, when it was simply encoding its own behavior. By regressing out each bat’s own behavior, this approach eliminates such potential spurious correlation between one bat’s activity and another bat’s behavior caused by coordinated behaviors between the bats. Figure 6A shows an example of one bat’s activity encoding the other bat’s behavior independently of its own behavior, and the strength of the encoding was quantified using the ROC curve in Figure 6B. Across simulations, we found that the encoding of the other bat’s behavior was significantly stronger for the onechamber model than the twochambers model (Figure 6C), consistent with the presence of functional acrossbrain coupling in the onechamber model. Applying the same approach to the data, we saw similar neural encoding of the other bat’s behavior (example in Figure 6D, corresponding ROC curve in Figure 6E). Across all data sessions, the encoding was stronger on onechamber sessions than twochambers sessions, confirming the model predictions (Figure 6F).
In summary, we tested and confirmed three disparate predictions of the model that all stem from the model mechanism of functional acrossbrain coupling. The fact that three seemingly unrelated experimental observations can all be explained by our model underscores its explanatory power.
Modeling interbrain neural relationship during group social interactions
So far, our analysis has been exclusively focused on social interactions between two individuals. However, social interactions also occur among larger groups in many species, including the groupliving Egyptian fruit bats (HerzigStraschil and Robinson, 1978; Kwiecinski and Griffiths, 1999). Although the interbrain relationship between two interacting individuals have been studied extensively, the interbrain relationship among larger social groups have received far less attention (but see Dikker et al., 2017; Rose et al., 2021). While there are many open questions regarding the group interbrain relationship, in the context of the current study, we ask two questions that are natural extensions of our twobat findings. First, can the opposite feedback mechanism be generalized to the larger group context? Second, what would such a mechanism predict about the group interbrain relationship? To answer these questions, we extended our model to larger groups, which offered direct predictions regarding the group interbrain relationship. We then tested these predictions using social group data in bats.
To extend the twobat model to interactions among a group of $n$ bats, we used the same equation shown in Figure 3A, except that the activity $\mathit{a}$ and behavioral modulation $\mathit{b}$ are now $n$dimensional rather than 2dimensional, and the functional coupling matrix $\mathit{C}$ is similarly $n\times n$. As in the twobat model, the activity of each bat in the $n$bat model is influenced by the same functional selfcoupling ${\mathit{C}}_{S}$ , and the activity of each pair of bats are functionally coupled with the same positive coupling $\mathit{C}}_{I$.
The mechanism that governs the twobat model naturally generalizes to the $n$bat model (see Materials and methods for details). For the $n$bat model, the mean component corresponds to the mean activity across all bats. On the other hand, a difference subspace takes the place of the difference component of the twobat model. In $n$bat activity space, the difference subspace is the $\left(n1\right)$dimensional subspace orthogonal to the direction of the mean component (Figure 7B). This subspace contains all interbrain activity patterns that correspond to activity differences across any of the brains. Similar to the twobat model, functional acrossbrain coupling acts as positive and negative feedback to the $n$bat mean component and the difference subspace, respectively, amplifying and slowing down the $n$bat mean component while suppressing and speeding up activity patterns in the difference subspace.
The opposite feedback mechanism in the $n$bat model thus gives rise to predictions that can be tested experimentally. Specifically, one prediction is that activity patterns in the difference subspace would have smaller magnitudes than the $n$bat mean component: as we show using simulations of a 4bat model, the average variance of the difference subspace is smaller than the variance of the 4bat mean component (Figure 7C). Another prediction is that activity patterns in the difference subspace would have faster timescales than the $n$bat mean component: 4bat simulations show that the difference subspace has a higher average power spectral centroid than the 4bat mean component (Figure 7D). Finally, similar to the twobat model (see Figure 5A–B), correlation between activity variables depends on their functional coupling. In particular, the positive functional acrossbrain coupling between the activity of different bats would give rise to positive interbrain correlations for each pair of bats within the group; on the other hand, the mean component is not functionally coupled with activity patterns in the difference subspace, so they would therefore be uncorrelated. This prediction is illustrated in Figure 7E using simulations of a 4bat model.
Next, we tested these predictions using experimental data collected from a group of socially interacting bats. In this data set, neural activity was simultaneously and wirelessly recorded from the frontal cortices of four bats while they interacted in a chamber, under the same condition as the onechamber sessions described above (Rose et al., 2021). As with the twobat experiment, we focused on LFP power in the 30–150 Hz band. Applying the same analyses used for the 4bat model simulations in Figure 7C–E, we found that the data qualitatively confirmed all model predictions: activity patterns in the difference subspace had smaller magnitudes (Figure 7F) and faster timescales (Figure 7G) than the 4bat mean component, there was positive interbrain correlation between bats in the group (Figure 7H), and there was no correlation between the 4bat mean component and activity patterns in the difference subspace (Figure 7H). Note that the 4bat experimental data and the model predictions differ quantitatively, because the 4bat model did not include realistic behavioral modeling as did the twobat model (Materials and methods).
In conclusion, we modeled the group interbrain relationship and experimentally confirmed the model predictions. Through this, we found that the opposite feedback mechanism is not restricted to pair social interactions, but naturally generalizes to larger social groups.
Discussion
Here we studied the interbrain relationship between socially interacting bats, by analyzing the mean and difference between the bats’ neural activity. We found that the difference in activity has smaller magnitudes and faster timescales than the mean, which is a robust neural signature of social interactions. We reproduced this finding using a simple model, which suggests a specific computational mechanism shaping the interbrain neural relationship: functional coupling across brains acting as opposite feedback to the difference and mean activity components. This mechanism gave rise to a number of predictions, which we then confirmed using experimental data from bats engaging in paired and group social interactions. Furthermore, this mechanism has broad explanatory power, and can account for a range of previously observed features of the interbrain neural relationship, as we discuss below.
First, it is well known that neural activity is correlated between the brains of socially interacting individuals, in mice, bats, and humans (e.g. Dikker et al., 2014; Kingsbury et al., 2019; Kinreich et al., 2017; Levy et al., 2017; Montague et al., 2002; Piazza et al., 2020; Spiegelhalder et al., 2014; Stolk et al., 2014; Zadbood et al., 2017; Zhang and Yartsev, 2019). Through negative feedback to the difference component, functional acrossbrain coupling suppresses the interbrain difference and keeps it small. At the same time, positive feedback to the mean component amplifies the common activity covariation between brains, so that the neural activity in each is dominated by their common activity pattern. Together, these would result in the interbrain correlation observed during social interactions in diverse species.
Second, in a shared social environment, neural correlation across brains is above and beyond what would be expected from behavioral coordination between individuals (Kingsbury et al., 2019; Piazza et al., 2020; Zhang and Yartsev, 2019). This can be understood because, through opposite feedback to the difference and mean components, functional acrossbrain coupling contributes to correlation independently from behavioral coordination. As a consequence of this, even in the absence of behavioral coordination, the opposite feedback results in correlation across brains in the model (Figure 4), consistent with previous experimental observations (Kingsbury et al., 2019; Piazza et al., 2020; Zhang and Yartsev, 2019).
Third, the neural correlation between socially interacting individuals varies as a function of timescale: correlation across brains is higher for activity at slower timescales, and lower for activity at faster timescales (Kingsbury et al., 2019; Zhang and Yartsev, 2019). This would be a natural consequence of the opposite feedback mechanism. Positive feedback to the mean amplifies it and slows it down, while negative feedback to the difference suppresses it and speeds it up. Thus, at slower timescales, the activity of the two brains is dominated by the mean component, that is, the activity of the two brains are very similar at slower timescales and thus have higher correlations. On the other hand, at faster timescales, the increasing presence of the difference component lowers the correlation across brains.
In summary, our model offers a particularly simple explanation for diverse aspects of the interbrain neural relationship during both pair and group social interactions, including: (1) neural correlation across brains and the relative magnitudes of interbrain mean and difference components (Figure 3E; Figure 7F and H); (2) neural correlation beyond behavioral coordination; (3) neural correlation as a function of timescale; (4) the relative timescales of interbrain mean and difference components (Figure 3F; Figure 7G); (5) correlation between activity variables under rotated activity bases (Figure 5B–C); (6) the relationship between the magnitudes and timescales of activity components across sessions (Figure 5G–H); and (7) modulation of one bat’s neural activity by the behavior of another bat (Figure 6). Importantly, there is no reason a priori to suspect that these disparate observations are all related; only in light of the model is it apparent that they are all manifestations of a single underlying computational mechanism. Collectively, these results thus suggest opposite feedback to the difference and mean as a unifying computational mechanism governing the interbrain relationship during social interactions. This points to a major avenue of future research: elucidating the biological implementation of this computational mechanism. Such research will likely involve studying the distributed neural circuits mediating sensorimotor interactions and attentional processes during social interactions, to understand how they might functionally act as opposite feedback to the interbrain difference and mean components. To do so, a relevant experiment is to record from multiple brain regions and examine how signals related to social interactions are transformed across regions. This will also reveal the extent to which interbrain correlation varies across regions, allowing comparison to the human literature, where variations across brain structures have long been studied using EEG and fMRI.
Another important future direction is comparative studies of the interbrain relationship to elucidate similarity and differences across species. So far, this has been studied in three nonhuman species: bats, mice (Kingsbury et al., 2019), and birds (Hoffmann et al., 2019). While bats and mice showed qualitatively similar interbrain relationships, birds were different, showing anticorrelated neural activity across brains during a duetting behavior (Hoffmann et al., 2019). It would be interesting to study the interbrain relationship in species with different levels of sociality and different types of social behaviors, for exmple, solitary animals that only socially interact for mating, or animals that live in large groups. Such comparative studies, with careful behavioral quantification, would identify neural mechanisms that underlie social interactions in general versus speciesspecific social behaviors. Moreover, the analytic and modeling methods used here can be easily applied to other species, enabling comparison across species under the same framework. In particular, studying the timescales of interbrain activity patterns during different social behaviors in different species could lead to a more general understanding of functional acrossbrain coupling.
Further future directions are motivated by the limitations of our model. For simplicity, our model took functional acrossbrain coupling to be static, that is, taking the same functional form and having the same strength independent of the bats’ behaviors. In reality, the coupling is mediated through behavioral interaction between the bats. Thus, it is likely dynamic, taking different forms and varying in strength depending on the bats’ dynamic behaviors. For example, when one bat is actively grooming another passively resting bat, functional coupling might be asymmetric between them, stronger in one direction than the other. Conversely, when two bats are fighting with each other, their functional coupling might be more symmetric. Furthermore, the timescales of functional coupling could also vary: different behaviors might couple neural activity with different time courses and delays. Future study focusing on specific behaviors (e.g. grooming, fighting, etc.), with finetimescale behavioral tracking (Mathis et al., 2018; Pereira et al., 2019), will be needed to identify the detailed behaviorspecific dynamics of functional acrossbrain coupling.
Another limitation of our model is the ‘openloop’ nature of the relationship between behavior and neural activity. Specifically, we modeled neural activity as being modulated by behavior, but behavior was modeled using a Markov chain that is independent from the neural activity. In reality, neural activity and behavior form a closedloop, with different social behaviors being controlled by the neural activity of specific neural populations in specific brain regions. Thus, an important future direction is to close the loop by incorporating neural control of social behaviors into models of the interbrain relationship in bats. This will require future experimental studies to identify which frontal cortical regions and populations in bats are necessary or sufficient to control social behaviors, as well as the detailed causal relationship from neural activity to social behavior. Furthermore, as social interactions can occur at multiple timescales, it will be interesting to investigate how these are controlled by neural activity at different timescales, and how those timescales are shaped by functional acrossbrain coupling. In summary, such a closedloop model will shed light on how interbrain activity patterns and dynamic social interactions coevolve and feedback onto each other.
Materials and methods
Request for resources
Request a detailed protocolFurther information and requests for resources should be directed to Michael M. Yartsev (myartsev@berkeley.edu).
Experimental and data processing methods
The data analyzed in this study have been previously published in Zhang and Yartsev, 2019 and Rose et al., 2021. The details of experimental methods and data processing methods for the twobat social interaction experiments were published in Zhang and Yartsev, 2019; the details for the group social interaction experiments were published in Rose et al., 2021; therefore, here we present a brief summary of the published methods. All data processing was performed using MATLAB (MathWorks).
Experimental subjects
Request a detailed protocolData was collected from eight adult male Egyptian fruit bats (Rousettus aegyptiacus). All procedures were approved by the Animal Care and Use Committee of UC Berkeley.
Experimental setup
Request a detailed protocolExperiments were conducted inside 40.6 × 33.7 × 52.1 cm (length × width × height) cages, placed in dark chambers. Highspeed infrared video cameras (Flea3 or Chameleon3, FLIR) were used to record videos of the bats, and ultrasonic microphones (USG Electret Ultrasound Microphone, Avisoft Bioacoustics) were used to record audio.
The experiments included onechamber sessions and twochambers sessions. In the onechamber sessions, bats behaved freely and interacted with each other inside a chamber. In the twochambers sessions, bats behaved freely in separate, identical chambers. There were three types of twochambers sessions: (1) two bats each freely behaving in isolation; (2) two bats each freely behaving in the presence of identical auditory stimuli (playback of bat calls); (3) two bats each freely behaving and interacting with a different partner in separate chambers.
Behavior annotation
Request a detailed protocolFor the twobat experiments, the behaviors of the bats were manually annotated by experienced observers. The annotated behaviors and their definitions are as follows.
Resting. A bat hanging by its feet, with its head and body still.
Active nonsocial. A bat engaging in any kind of active behavior that does not involve social interaction, including: the bat hanging by its feet or feet and thumbs, and moving its head or body; the bat climbing or crawling around; the bat shaking its body; the bat jumping or flying off from the roof of the cage.
Selfgrooming. A bat either licking or scratching itself.
Social grooming. A bat either licking or scratching another bat.
Probing. A bat poking its snout at the head or body of another bat.
Fighting. A bat moving its wings or thumbs to quickly hit another bat, or biting another bat.
Mating. A male bat inserting or attempting to insert its penis into a female bat’s vagina.
Wing covering. A bat struggling with another bat in order to cover the other bat’s body with its opened wings.
Reaching. A bat attempting to reach over the body or wings of another bat with its head or thumbs.
Blocking. A bat using its wings to actively block another bat from accessing a location.
Other interactions. Any social interaction other than the ones already defined.
Surgery
Request a detailed protocolAnesthesia and surgery were performed as described previously in Yartsev and Ulanovsky, 2013 to implant the bats with fourtetrode microdrives (Harlan 4 Drive, Neuralynx). To target the frontal cortex, the center of craniotomy was positioned at 1.7 mm lateral to the midline and 12.19 mm anterior to the transverse sinus located between cortex and cerebellum.
Electrophysiological recording
Request a detailed protocolElectrophysiology was performed using a wireless neural data logging system (Neurolog16 and MouseLog16, Deuteron Technologies). Tetrodes were advanced ventrally every one to two days (mostly by 20–160 µm), to record activity at different sites.
Histology
Request a detailed protocolTo determine the neural recording sites, we performed histology as previously described (Yartsev and Ulanovsky, 2013), using Nissl staining.
Preprocessing of electrophysiological data
Request a detailed protocolTo obtain LFP, we lowpass filtered the raw voltage traces (cutoff frequency: 200 Hz), and then downsampled it to 496.6 Hz (the twobat experiments) or 520.8 Hz (the fourbat experiments).
We detected spikes from bandpass filtered (600–6000 Hz) voltage traces using threshold crossing. Spike sorting was done automatically using SNAP Sorter (Neuralynx), then manually checked using SpikeSort3D (Neuralynx). For each tetrode on each session, all spikes not assigned to single units were grouped into a multiunit. All units with firing rate below 2 Hz were excluded from further analysis.
Calculation of LFP spectrograms
Request a detailed protocolTo calculate LFP spectrograms, power spectra were calculated using the multitaper method for 5 s sliding windows of the LFP trace, with 2.5 s overlap between consecutive windows. To analyze different frequencies of the LFP on equal footing, we separately peaknormalized the power at each frequency for each spectrogram.
Calculation of firing rates
Request a detailed protocolFiring rates for single units and multiunits were computed in 5 s bins with 2.5 s overlap between consecutive bins.
Dimensionality reduction of LFP
Request a detailed protocolPreviously, we found that the spectrograms of bat frontal cortical LFP can be reduced to two signals, power in the 1–29 Hz band and the 30–150 Hz band (Zhang and Yartsev, 2019). To analyze these signals, for a given normalized LFP spectrogram (with power peaknormalized for each frequency), we averaged the normalized power values across 1–29 Hz and across 30–150 Hz at each time bin. This is the ‘mean normalized LFP power’ used for all analyses in this paper involving LFP power (e.g. Figure 1C).
Analysis and modeling
All analysis and modeling were performed using MATLAB (MathWorks).
Statistical tests
Request a detailed protocolThe statistical tests used are stated in the figure legends. A significance level of 0.05 was used for all tests. Tests were twotailed unless otherwise indicated.
Interbrain difference and mean components: definition and quantification
Request a detailed protocolThe neural activity of two bats can be represented as a twodimensional vector $\left(\begin{array}{c}{\mathit{a}}_{1}\left(t\right)\\ {\mathit{a}}_{2}\left(t\right)\end{array}\right)$, where ${\mathit{a}}_{1}\left(t\right)$ and ${\mathit{a}}_{2}\left(t\right)$ are the neural activity (normalized LFP power, multiunit activity, or single unit activity) of bat 1 and bat 2 at time $t$, respectively. Through a change of basis, the same activity can be represented under another orthogonal basis as the mean and difference between the two brains: $\left(\begin{array}{c}{\mathit{a}}_{M}\left(t\right)\\ {\mathit{a}}_{D}\left(t\right)\end{array}\right)=\frac{1}{2}\left(\begin{array}{c}{\mathit{a}}_{1}\left(t\right)+{\mathit{a}}_{2}\left(t\right)\\ {\mathit{a}}_{1}\left(t\right){\mathit{a}}_{2}\left(t\right)\end{array}\right)$, where ${\mathit{a}}_{M}\left(t\right)$ is the mean component of the activity, and ${\mathit{a}}_{D}\left(t\right)$ is the difference component. Here, the difference component is defined as $\frac{1}{2}\left[{\mathit{a}}_{1}\left(t\right){\mathit{a}}_{2}\left(t\right)\right]$ rather than ${\mathit{a}}_{1}\left(t\right){\mathit{a}}_{2}\left(t\right)$ so as to have the same scale as the mean component. ${\mathit{a}}_{M}\left(t\right)$ represents the common activity pattern shared between the brains, and ${\mathit{a}}_{D}\left(t\right)$ represents the activity difference between the brains.
To quantify the magnitude of an activity component on a given session, we computed the variance of its time series over the session. To quantify the timescales of an activity component on a given session, we computed its power spectral centroid as follows. Given the time series of an activity component over the session, we subtracted from it its average over time, multiplied it with a Hamming window, and then computed the periodogram estimate of its power spectrum. The power spectral centroid is then computed as a weighted average of frequency, with the power at each frequency as its weight. A higher power spectral centroid means that the activity has more power at higher frequencies, so that the time series varied at faster timescales. Note that the power spectral centroid was always calculated from time series of mean normalized LFP power or firing rate, and not from time series of LFP itself.
Surrogate data and the relationship between interbrain correlation and mean and difference components
Request a detailed protocolIn our data from onechamber sessions, the mean component had large magnitudes and slow timescales, and the difference component had small magnitudes and fast timescales. The relative magnitudes of the two components are expected, given that activity from the two bats showed similar patterns over time (Figure 1C), with high positive interbrain correlations (Zhang and Yartsev, 2019); as explained below, a positive correlation mathematically implies larger variance for the mean component compared to the difference component. What about their timescales? Are the observed relative timescales of the two components necessary mathematical consequences of high interbrain correlations, large mean components, and small difference components?
To answer this, we examine the relationship between interbrain correlation and the mean and difference components. Let’s use $\mathit{A}}_{1$, a column vector, to denote the activity of bat 1 during a session, and similarly, $\mathit{A}}_{2$ for the activity of bat 2. $\mathit{A}}_{1$ and $\mathit{A}}_{2$ are $N$dimensional vectors, where $N$ is the number of time points in the session. The activity of the mean component is ${\mathit{A}}_{M}=\left({\mathit{A}}_{1}+{\mathit{A}}_{2}\right)/2$, and the activity of the difference component is ${\mathit{A}}_{D}=\left({\mathit{A}}_{1}{\mathit{A}}_{2}\right)/2$. We use $\hat{\mathit{A}}}_{1$ to denote the meansubtracted version of $\mathit{A}}_{1$ (here ‘meansubtracted’ refers to subtracting the average across the elements of a vector, that is, across time, not to be confused with the average across bats, as in the ‘mean component’), and similarly, $\hat{\mathit{A}}}_{2$, $\hat{\mathit{A}}}_{M$, and $\hat{\mathit{A}}}_{D$ for meansubtracted versions of the respective vectors. Then, the Pearson correlation coefficient between the activity of two brains is:
Thus, the correlation depends on the quantities $\hat{\mathit{A}}}_{M}^{T}{\hat{\mathit{A}}}_{M$, $\hat{\mathit{A}}}_{D}^{T}{\hat{\mathit{A}}}_{D$, and $\left({\hat{\mathit{A}}}_{M}^{T}{\hat{\mathit{A}}}_{D}\right)}^{2$. Note that $\hat{\mathit{A}}}_{M}^{T}{\hat{\mathit{A}}}_{M$ is $N$ times the variance of the mean component, and $\hat{\mathit{A}}}_{D}^{T}{\hat{\mathit{A}}}_{D$ is $N$ times the variance of the difference component. A positive correlation requires the numerator in equation (1) to be positive, in other words, the mean component having larger variance than the difference component, as stated above.
Equation (1) also shows that having a given combination of correlation, mean component variance, and difference component variance does not place constraints on the timescales of the mean and difference components. Specifically, it does not constrain the difference component to be faster than the mean. To explicitly demonstrate this, we generated surrogate data with identical interbrain correlation, mean component variance, and difference component variance as actual data, but having a slower difference component than the mean component (Figure 2).
The surrogate data in Figure 2A–C were generated from the actual 30–150 Hz LFP power data of the example session shown in Figure 1C, by keeping the mean component of the actual data, and replacing the difference component with a surrogate, using the following procedure. Below we denote the actual data by $\mathit{A}}_{D$, $\hat{\mathit{A}}}_{D$, etc. as above, and denote their surrogate counterparts by $\mathit{S}}_{D$, $\hat{\mathit{S}}}_{D$, etc. We generated a random $N\times 1$ activity vector by picking each element independently from the uniform distribution between 0 and 1, smoothed it with a 1000second moving average filter, and subtracted from it the average across its elements. Let’s call the resulting vector $\hat{\mathit{R}}$. Due to the smoothing, $\hat{\mathit{R}}$ varied on slow timescales. We then found the component of $\hat{\mathit{R}}$ that is orthogonal to $\hat{\mathit{A}}}_{M}:{\hat{\mathit{R}}}_{O}=\hat{\mathit{R}}[({\hat{\mathit{R}}}^{T}{\hat{\mathit{A}}}_{M})/({\hat{\mathit{A}}}_{M}^{T}{\hat{\mathit{A}}}_{M})]{\hat{\mathit{A}}}_{M$. Next, we constructed a vector $\hat{\mathit{S}}}_{D}^{\prime$ as a linear combination of $\hat{\mathit{A}}}_{M$ and ${\hat{\mathit{R}}}_{O}:{\hat{\mathit{S}}}_{D}^{\prime}=$ $\mathrm{c}\mathrm{o}\mathrm{t}[\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{c}\mathrm{o}\mathrm{s}({\hat{\mathit{A}}}_{M}^{T}{\hat{\mathit{A}}}_{D}/\sqrt{{\hat{\mathit{A}}}_{M}^{T}{\hat{\mathit{A}}}_{M}{\hat{\mathit{A}}}_{D}^{T}{\hat{\mathit{A}}}_{D}})]{\hat{\mathit{A}}}_{M}/\sqrt{{\hat{\mathit{A}}}_{M}^{T}{\hat{\mathit{A}}}_{M}}$ $+{\hat{\mathit{R}}}_{O}/\sqrt{{\hat{\mathit{R}}}_{O}^{T}{\hat{\mathit{R}}}_{O}}$. The surrogate for $\hat{\mathit{A}}}_{D$ is the vector $\hat{\mathit{S}}}_{D$, obtained by scaling $\hat{\mathit{S}}}_{D}^{\prime$ to have the same vector norm as $\hat{\mathit{A}}}_{D}:{\hat{\mathit{S}}}_{D}=(\sqrt{{\hat{\mathit{A}}}_{D}^{T}{\hat{\mathit{A}}}_{D}}/\sqrt{{\hat{\mathit{S}}}_{D}^{\prime T}{\hat{\mathit{S}}}_{D}^{\prime}}){\hat{\mathit{S}}}_{D}^{\prime$. The surrogate difference component is then $\mathit{S}}_{D}={\hat{\mathit{S}}}_{D}+{\overline{\mathit{A}}}_{D$, where every element of the vector $\overline{\mathit{A}}}_{D$ is equal to the average across the elements of $\mathit{A}}_{D$. The surrogate activity of bat 1 and bat 2 are, respectively, $\mathit{S}}_{1}={\mathit{A}}_{M}+{\mathit{S}}_{D$ and $\mathit{S}}_{2}={\mathit{A}}_{M}{\mathit{S}}_{D$. This procedure ensures that $\hat{\mathit{S}}}_{D}^{T}{\hat{\mathit{S}}}_{D}={\hat{\mathit{A}}}_{D}^{T}{\hat{\mathit{A}}}_{D$ and $\left({\hat{\mathit{A}}}_{M}^{T}{\hat{\mathit{S}}}_{D}\right)}^{2}={\left({\hat{\mathit{A}}}_{M}^{T}{\hat{\mathit{A}}}_{D}\right)}^{2$, so that the surrogate data had identical interbrain correlation, mean component variance, and difference component variance as the actual data from which it was generated. Note that here we chose to leave the mean component $\mathit{A}}_{M$ from the actual data unchanged in generating the surrogate data, so that the surrogate $\mathit{S}}_{1$ and $\mathit{S}}_{2$ show qualitatively similar behavioral modulation over time as the actual data $\mathit{A}}_{1$ and $\mathit{A}}_{2$, but it is also possible to replace both $\mathit{A}}_{M$ and $\mathit{A}}_{D$ with surrogate counterparts.
For Figure 2D–E, we repeated the above procedure for each of the onechamber sessions from the actual data set, using the actual 30–150 Hz LFP power data.
Modeling: procedure
Request a detailed protocolOur goal of modeling was to reproduce the observed relationship between difference and mean components using a model that is simple and parsimonious, in order to unambiguously identify the underlying computational mechanisms. We modeled the time evolution of the neural activity of two bats using the following differential equation:
Here, $\tau $ is a time constant; $t$ is time; $\mathit{a}\left(t\right)=\left(\begin{array}{c}{\mathit{a}}_{1}\left(t\right)\\ {\mathit{a}}_{2}\left(t\right)\end{array}\right)$ is the activity of bat 1 and bat 2 at time $t$; $\mathit{C}=\left(\begin{array}{cc}{\mathit{C}}_{S}& {\mathit{C}}_{I}\\ {\mathit{C}}_{I}& {\mathit{C}}_{S}\end{array}\right)$ is the functional coupling matrix, where $\mathit{C}}_{S$ is the strength of functional selfcoupling and $\mathit{C}}_{I$ is the strength of functional acrossbrain coupling; $\mathit{b}\left(t\right)=\left(\begin{array}{c}{\mathit{b}}_{1}\left(t\right)\\ {\mathit{b}}_{2}\left(t\right)\end{array}\right)$ is the strengths of behavioral modulation, where ${\mathit{b}}_{1}\left(t\right)$ is the modulation of the neural activity of bat 1 by the behavior of bat 1 at time $t$, and similarly for ${\mathit{b}}_{2}\left(t\right)$ .
The functional acrossbrain coupling term, $\mathit{C}}_{I$, represents the indirect influence (as opposed to direct influence from an actual neural connection) one bat’s neural activity has on the other bat’s neural activity. For onechamber sessions, ${\mathit{C}}_{I}>0$, which models positive functional coupling when the two bats share a social environment (Hasson et al., 2012). For example, when bat 1’s neural activity (say, 30–150 Hz LFP power) increases due to its active movements (Gervasoni et al., 2004; McGinley et al., 2015), the movements create sensory inputs to bat 2, which can increase bat 2’s neural activity to the extent that bat 2 is paying attention (Chun et al., 2011; Driver, 2001; Fritz et al., 2007; Reynolds and Chelazzi, 2004). For twochambers sessions, ${\mathit{C}}_{I}=0$ since the two bats are separated. To ensure stability (so that neural activity do not go to infinity), ${\mathit{C}}_{S}$ must be negative and must have a larger absolute value than $\mathit{C}}_{I$ ; thus, $0\le {\mathit{C}}_{I}<{\mathit{C}}_{S}$ .
To generate $\mathit{b}\left(t\right)$ for a simulation, we first simulated the bats’ behaviors using a Markov chain. Each state of the Markov chain corresponds to the behaviors of the two bats at a given time: for example, bat 1 resting and bat 2 selfgrooming would be one state. For the Markov chain for onechamber sessions, the transition probability matrix, initial distribution, and state space were determined as follows. The transition probability from one state to another is taken to be the empirical frequency of that transition during all onechamber sessions. To calculate the empirical frequency, the behavior of each bat was sampled every 2.5 s, at the same time points as the neural activity (i.e. the center time point of each window used to calculate LFP power). A transition was counted for each consecutive pair of time points (in the rare instances where a bat engaged in multiple behaviors at the same time point, one transition was counted for each of the simultaneous behaviors). The initial distribution was taken to be the distribution of the behavioral state at the first time point of each onechamber session. For these calculations, the two bats were assumed to be symmetrical, in the following sense. Using ‘AB’ to denote the state of ‘bat 1 engaging in behavior A, bat 2 engaging in behavior B’, we consider the transitions ‘AB→CD’ and ‘BA→DC’ to have the same transition probability. When calculating the empirical transition frequencies, the count for ‘AB→CD’ and the count for ‘BA→DC’ were each taken to be the sum of the actual counts for ‘AB→CD’ and ‘BA→DC’. This applied also to states where both bats were engaging in the same behavior: for exmple, the count for the transition ‘AA→AA’ was taken to be double the actual count. This symmetry assumption was made for simplicity, and to allow behavioral data from different pairs of bats to be pooled together. Once the empirical transition frequencies were calculated, if there were less than 100 transitions from a given state, then that state was excluded from the state space. The same procedures were used to determine the transition probability matrix, initial distribution, and state space for the twochambers Markov chain. The transition probability matrices for onechamber and twochambers sessions are shown in Figure 3—figure supplement 2CD.
For each simulated session, we simulated the bats’ behaviors for 100 minutes using the Markov chain (with 2.5 s between steps, the behavioral sampling period used when calculating the empirical transition frequencies). The first state of each simulated session was drawn randomly from the appropriate initial distribution. From the simulated behaviors, $\mathit{b}\left(t\right)$ was determined at the discrete time points of the Markov chain steps (at integer multiples of 2.5 s). At the time point of a given Markov chain step, say ${t}_{1}$, $\mathit{b}\left({t}_{1}\right)$ was the sum of the noiseless deterministic behavioral modulation ${\mathit{b}}_{d}\left({t}_{1}\right)$ and the behavioral modulation noise ${\mathit{b}}_{n}\left({t}_{1}\right)$. ${\mathit{b}}_{d}\left({t}_{1}\right)$ depended on the behaviors of the two bats at time ${t}_{1}$. For example, if bat 1 was resting and bat 2 was engaging in selfgrooming at time ${t}_{1}$, then $\mathit{b}}_{d}\left({t}_{1}\right)=\left(\begin{array}{c}{b}_{resting}\\ {b}_{selfgrooming}\end{array}\right)+{b}_{constant$, where ${b}_{resting}$ and ${b}_{selfgrooming}$ are parameters specifying the level of behavioral modulation associated with the two behaviors, and ${b}_{constant}$ is a constant offset that differs between the onechamber and twochambers models, reflecting the effects of the general level of arousal that differs between the two conditions. See section ‘Modeling: parameters’ below for the parameter values for all the behaviors and ${b}_{constant}$. For the behavioral modulation noise, the two elements (for the two bats) of ${\mathit{b}}_{n}\left({t}_{1}\right)$ was each drawn independently from a Gaussian distribution with zero mean and standard deviation ${\sigma}_{n}$. After determining $\mathit{b}\left(t\right)={\mathit{b}}_{d}\left(t\right)+{\mathit{b}}_{n}\left(t\right)$ at the discrete time points of the Markov chain steps, $\mathit{b}\left(t\right)$ at time points inbetween were linearly interpolated.
Having simulated the behaviors and generated $\mathit{b}\left(t\right)$ for a 100 minute session, we set the initial condition $\mathit{a}\left(0\right)$ to be the fixed point under the initial noiseless behavioral modulation ${\mathit{b}}_{d}\left(0\right)$: $\mathit{a}\left(0\right)={\mathit{C}}^{1}{\mathit{b}}_{d}\left(0\right)$. Then, equation (2) was numerically integrated (ode45 function in MATLAB) to simulate neural activity $\mathit{a}\left(t\right)$. The only differences between simulations of onechamber and twochambers sessions were the value of $\mathit{C}}_{I$ and the behavior transition probability matrix. For analyses and figures involving $\mathit{a}\left(t\right)$, $\mathit{b}\left(t\right)$, ${\mathit{b}}_{d}\left(t\right)$, and ${\mathit{b}}_{n}\left(t\right)$ (Figures 3—7 and their figure supplements), we used their values taken at discrete time points corresponding to the Markov chain steps (i.e. same sampling period as the data).
Modeling: comparing behavior models
Request a detailed protocolAs described in the previous section, we modeled the behavior of the bats using a Markov chain. To justify this choice, we tested the Markov assumption by comparing three behavior models, including the Markov chain model we used. The first model is the independent model, where the behavioral state at a given time point is independent from the state at other time points. The second model is the firstorder dependency model, where the behavioral state at a given time point depends on the state at the previous time point only; this was implemented as a part of our main model in the previous section, and it corresponds to the Markov assumption. The third model is the secondorder dependency model, where the behavioral state at a given time point depends on the states at the two previous time points. Models with longer timedependencies (≥ 3) were not tested because the number of parameters grows exponentially with model order and cannot be fitted with our dataset.
The three models were each fitted separately for the onechamber and twochambers sessions. Each type of sessions was split into a training set and a test set of sessions (80% of the sessions for training and 20% for test). The models were fit by setting the probabilities or conditional probabilities of behaviors to the empirically observed frequencies (Laplace smoothing was used to avoid assigning zero probability to unobserved events). Then, the loglikelihood of the test set under each model was calculated. This procedure was repeated for 100 random splits of the data into training and test sets, and the loglikelihoods are shown in Figure 3—figure supplement 2A, B.
As shown in Figure 3—figure supplement 2A, B, the firstorder model had the highest likelihood on average. This does not necessarily prove that bat behavior obeys the Markov assumption (e.g. given more data, it is possible that better secondorder and higher order models could be fit). But it does mean that, given the amount of data we have, the best behavior model that we can fit is the firstorder Markov chain, supporting its use in our main model described in the previous section.
Modeling: comparing actual and simulated neural activity as a function of time
Request a detailed protocolTo compare neural activity as a function of time between the data and the model, we simulated our model using experimentally observed behaviors. In other words, rather than generating $\mathit{b}\left(t\right)$ from behaviors simulated using a Markov chain, we generated $\mathit{b}\left(t\right)$ based on the actual behaviors observed on a given experimental session. Neural activity was then simulated and compared to the actual neural activity on that experimental session. To quantify how well the model reproduces neural activity over time, we calculated the correlation between simulated and actual neural activity for each session and each bat: the average correlation over all sessions and bats is 0.72 (standard deviation 0.10). An example session is shown in Figure 3—figure supplement 1.
Modeling: magnitudes and timescales of mean and difference components
Request a detailed protocolHow is the model able to reproduce the set of experimental observations relating the magnitudes and timescales of the mean and difference components (Figure 3E–F)? To understand the mechanisms in the model, we examine equation (2). In general, in equations of this form, the activity variables are coupled: for example, in the onechamber model, the activity of bat 1 ${\mathit{a}}_{1}\left(t\right)$ influences the activity of bat 2 ${\mathit{a}}_{2}\left(t\right)$, which in turn feeds back on ${\mathit{a}}_{1}\left(t\right)$. It is easier to understand how the activity evolves if we express the equation in a basis in which the activity variables uncouple. This basis consists of the eigenvectors of the functional coupling matrix $\mathit{C}$. The eigenvectors of $\mathit{C}$ are $\left(\begin{array}{c}1\\ 1\end{array}\right)$ and $\left(\begin{array}{c}1\\ 1\end{array}\right)$, with respective eigenvalues ${\mathit{C}}_{S}+{\mathit{C}}_{I}$ and ${\mathit{C}}_{S}{\mathit{C}}_{I}$. We define $\mathit{V}$ to be the matrix whose columns are the eigenvectors: $\mathit{V}=\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right)$. Expressed in the eigenvector basis, the activity variables become the mean activity across bats, and the difference in activity between bats: ${\mathit{V}}^{1}\left(\begin{array}{c}{\mathit{a}}_{1}\left(t\right)\\ {\mathit{a}}_{2}\left(t\right)\end{array}\right)=\frac{1}{2}\left(\begin{array}{c}{\mathit{a}}_{1}\left(t\right)+{\mathit{a}}_{2}\left(t\right)\\ {\mathit{a}}_{1}\left(t\right){\mathit{a}}_{2}\left(t\right)\end{array}\right)=\left(\begin{array}{c}{\mathit{a}}_{M}\left(t\right)\\ {\mathit{a}}_{D}\left(t\right)\end{array}\right)$. Thus, the uncoupled activity variables are precisely what we are interested in: the mean and difference components. Similarly, in the eigenvector basis, the behavioral modulation variables are the mean behavioral modulation and the difference in behavioral modulation: ${\mathit{V}}^{1}\left(\begin{array}{c}{\mathit{b}}_{1}\left(t\right)\\ {\mathit{b}}_{2}\left(t\right)\end{array}\right)=\frac{1}{2}\left(\begin{array}{c}{\mathit{b}}_{1}\left(t\right)+{\mathit{b}}_{2}\left(t\right)\\ {\mathit{b}}_{1}\left(t\right){\mathit{b}}_{2}\left(t\right)\end{array}\right)=\left(\begin{array}{c}{\mathit{b}}_{M}\left(t\right)\\ {\mathit{b}}_{D}\left(t\right)\end{array}\right)$. Thus, in the eigenvector basis, equation (2) is
In the onechamber model, ${\mathit{C}}_{I}>0$, so functional acrossbrain coupling acts as positive feedback to the mean component and negative feedback to the difference component. The effects of such opposite feedback to the two components can be seen more clearly by rewriting equation (3) as the following uncoupled differential equations:
Thus, at any given time $t$, ${\mathit{a}}_{M}\left(t\right)$ is exponentially approaching $\frac{{\mathit{b}}_{M}\left(t\right)}{{\mathit{C}}_{S}{\mathit{C}}_{I}}$ with effective time constant $\frac{\tau}{{\mathit{C}}_{S}{\mathit{C}}_{I}}$, and ${\mathit{a}}_{D}\left(t\right)$ is exponentially approaching $\frac{{\mathit{b}}_{D}\left(t\right)}{{\mathit{C}}_{S}+{\mathit{C}}_{I}}$ with effective time constant $\frac{\tau}{{\mathit{C}}_{S}+{\mathit{C}}_{I}}$. In the onechamber model, $0<{\mathit{C}}_{I}<{\mathit{C}}_{S}$, so $0<{\mathit{C}}_{S}{\mathit{C}}_{I}<{\mathit{C}}_{S}+{\mathit{C}}_{I}$. This means that in the onechamber model, relative to the twochambers model where ${\mathit{C}}_{I}=0$, the positive feedback provided by functional acrossbrain coupling amplifies the mean component and slows it down. On the other hand, the negative feedback suppresses the difference component and speeds it up. Thus, the model suggests this opposite feedback to be a potential computational mechanism that explains all of our observations relating the magnitudes and timescales of the mean and difference components.
The only differences between the onechamber and twochambers models were the value of $\mathit{C}}_{I$ and the behavioral transition probability matrix. Having understood the effects of $\mathit{C}}_{I$ , we now turn to the effects of the different transition probability matrices, which govern the Markov chain behavior models. For each Markov chain, we calculated its stationary distribution as the eigenvector of its transition probability matrix corresponding to an eigenvalue of 1, which was confirmed numerically in Figure 3—figure supplement 2E, F. Figure 3—figure supplement 2E, F shows that, for both the onechamber and twochambers models, the respective Markov chains approach stationary distributions in ~5 minutes, that is, the distribution of states no longer changes with time after ~5 min. The stationary distribution for the twochambers model shows that states where the two bats engage in the same behavior are about equally likely as states where they engage in different behaviors (Figure 3—figure supplement 2G), as expected given that the two bats behave independently in separate chambers. For the onechamber model, samebehavior states (probability 0.58) are more likely than differentbehavior states (Figure 3—figure supplement 2G). The noiseless behavioral modulation ${\mathit{b}}_{d}\left(t\right)$ depends on the behaviors of the two bats, while the behavioral modulation noise ${\mathit{b}}_{n}\left(t\right)$ does not. Through ${\mathit{b}}_{d}\left(t\right)$, samebehavior states contribute to behavioral modulation of the mean component, and the only noiseless behavioral modulation of the difference component comes from differentbehavior states. In the onechamber model, with increased coordinated behavioral modulation compared to the twochambers model, the mean noiseless behavioral modulation across bats has larger magnitudes and slower timescales than its difference between bats (blue lines in Figure 3—figure supplement 2H, I). On the other hand, because the behavioral modulation noise is independent across bats, it does not distinguish between the mean and difference components. Thus, adding the noise to the noiseless behavioral modulation decreases the differential modulation of the mean and difference components (red lines in Figure 3—figure supplement 2H, I). However, even if the noise strongly reduces differential behavioral modulation of the mean and difference components, the experimentally observed levels of relative magnitudes and timescales can still result from functional acrossbrain coupling (purple lines in Figure 3—figure supplement 2H, I).
Thus, our model suggests two mechanisms behind our observations on the mean and difference components: opposite feedback by functional acrossbrain coupling, and coordinated behavioral modulation. Due to these mechanisms, in a shared social environment, the activity of the two bats become dominated by their common activity pattern, and when the activity of the two bats diverge, they rapidly converge again. As a result, the activity of two socially interacting bats are highly correlated with each other. We found that interbrain correlation remained even when the bats were engaged in different behaviors (Figure 4A–B; see also Zhang and Yartsev, 2019), which suggests that coordinated behavioral modulation alone is not sufficient to explain the data, and that opposite feedback by functional acrossbrain coupling is needed. To explicitly test this, we performed an additional set of simulations of the onechamber model. Here, after we generated each instantiation of $\mathit{b}\left(t\right)$, we used it to simulate neural activity twice: once with functional acrossbrain coupling and once without. Note that the same $\mathit{b}\left(t\right)$ is used for both simulations, including the same ${\mathit{b}}_{d}\left(t\right)$ generated from the onechamber Markov chain and the same instantiation of randomly generated ${\mathit{b}}_{n}\left(t\right)$. Figure 4C, D and G shows that the model with functional acrossbrain coupling reproduced the experimental observation: interbrain correlation remained after removing periods of coordinated behaviors. On the other hand, without functional acrossbrain coupling, interbrain correlation disappeared after removing coordinated behaviors (Figure 4E–G). Thus, we conclude that, in the framework of our model, opposite feedback by functional acrossbrain coupling is necessary to reproduce the data.
Modeling: mean and difference components in a reduced model
Request a detailed protocolTo gain a more precise understanding how the relative variances and timescales of the mean and difference components depend on model parameters, we studied a reduced version of the model. In the reduced model, instead of generating $\mathit{b}\left(t\right)$ based on simulated behavior, ${\mathit{b}}_{M}\left(t\right)$ and ${\mathit{b}}_{D}\left(t\right)$ are simply noise with identical, flat power spectra. Such noise amounts to inputs with the same variance and same timescales to both the mean and difference components. Thus, the relative variances and timescales of the two components are solely determined by the coupling parameters $\mathit{C}}_{S$ and $\mathit{C}}_{I$. In the following, we derive expressions for the variance ratio and the power spectral centroid ratio of the mean and difference components, as functions of the coupling parameters.
Rewriting equation (3) and assuming $\tau =1$ for simplicity (as $\tau $ is a redundant parameter and can be absorbed into the other parameters), the reduced model can be written as
The solutions to these equations are:
In the above, ${f}_{M}\left(t\right)$ and ${f}_{D}\left(t\right)$ are the neural decay functions for the mean and difference components:
where $\lambda}_{M}={\mathit{C}}_{S}+{\mathit{C}}_{I$ and $\lambda}_{D}={\mathit{C}}_{S}{\mathit{C}}_{I$ are the eigenvalues of $\mathit{C}$.
The convolution theorem can be used to approximate the solutions (6) and (7) in the frequency domain:
Here, $T$ is the duration of a simulated session, and the tilde denotes a Fourier coefficient, for exmple, ${\stackrel{~}{f}}_{M,k}$ is the Fourier coefficient of ${f}_{M}\left(t\right)$ at frequency $k$: $\stackrel{\sim}{f}}_{M,k}=\frac{1}{\sqrt{T}}{\int}_{0}^{T}dt\text{}{f}_{M}\left(t\right){e}^{i\frac{2\pi kt}{T}$. The approximation of (10) and (11) assumes periodic convolutions and $\mathit{a}\left(0\right)=0$ in (6) and (7), but it is still accurate when $\mathit{a}\left(0\right)$ is smaller or comparable to the scale of $\mathit{b}\left(t\right)$ and $T$ is large compared to the effective time constants $\frac{1}{{\mathit{C}}_{S}{\mathit{C}}_{I}}$ and $\frac{1}{{\mathit{C}}_{S}+{\mathit{C}}_{I}}$.
We first consider the variance ratio. Using Parseval’s theorem and (10), the variance of the mean component can be approximated in the frequency domain as
In the reduced model, $\mathit{b}}_{M$ has a flat power spectrum: ${\left{\stackrel{\sim}{\mathit{b}}}_{M,k}\right}^{2}=p$ for all $k\ge 1$. Then, ${Var}_{M}\approx 2p{\sum}_{k=1}^{\infty}{\left{\stackrel{~}{f}}_{M,k}\right}^{2}$. Transforming back to the time domain, plugging (8) in, and evaluating the integrals, we have
The same calculation applied to the variance of the difference component shows that it is
Assuming ${\lambda}_{M}<0$ and ${\lambda}_{D}<0$ (required for stability), for large $T$, (12) and (13) can be approximated as
Thus, the variance ratio (variance of the mean divided by variance of the difference) is
Next, we consider the power spectral centroid ratio. Using (10), the power spectral centroid of the mean component can be approximated as
The Fourier coefficients of ${f}_{M}\left(t\right)$ are ${\stackrel{\sim}{f}}_{M,k}=\frac{1}{\sqrt{T}}{\int}_{0}^{T}dt\phantom{\rule{thinmathspace}{0ex}}{e}^{{\lambda}_{M}t}{e}^{i\frac{2\pi kt}{T}}=\frac{1}{\sqrt{T}}\frac{1}{{\lambda}_{M}i\frac{2\pi k}{T}}\left({e}^{{\lambda}_{M}T}1\right)$. Plugging this into (16), we have:
where ${\alpha}_{M}=\frac{{T}^{2}{\lambda}_{M}^{2}}{4{\pi}^{2}}$. The series ${\sum}_{k=1}^{\infty}\frac{1}{\frac{{\alpha}_{M}}{k}+k}$ can be shown to be divergent by comparison to $\frac{1}{{\alpha}_{M}+1}{\sum}_{k=1}^{\infty}\frac{1}{k}$, the harmonic series scaled by $\frac{1}{{\alpha}_{M}+1}$: $\frac{1}{\frac{{\alpha}_{M}}{k}+k}\ge \frac{1}{{\alpha}_{M}k+k}$ for all $k\ge 1$.
The same calculation applied to the power spectral centroid of the difference component shows that it is
where $\alpha}_{D}=\frac{{T}^{2}{\lambda}_{D}^{2}}{4{\pi}^{2}$.
While the power spectral centroids involve divergent series, we can define the power spectral centroid ratio (centroid of the mean divided by centroid of the difference) using a limit of the ratio of partial sums:
For simplicity, we rewrite the limit in (17) as $\underset{n\to \infty}{\mathit{lim}}\frac{{\beta}_{n}^{M}}{{\beta}_{n}^{D}}$ where ${\beta}_{n}^{M}={\sum}_{k=1}^{n}\frac{1}{\frac{{\alpha}_{M}}{k}+k}$ and ${\beta}_{n}^{D}={\sum}_{k=1}^{n}\frac{1}{\frac{{\alpha}_{D}}{k}+k}$. Note that ${\beta}_{n}^{D}$ is strictly increasing and $\underset{n\to \infty}{\mathit{li}m}{\beta}_{n}^{D}=\infty $, so we can evaluate the limit in (17) using the Stolz–Cesàro theorem:
Thus, assuming ${\lambda}_{M}<0$ and ${\lambda}_{D}<0$, for large $T$, (17) can be approximated as
The results $r}_{Var}\approx \frac{{\mathit{C}}_{S}+{\mathit{C}}_{I}}{{\mathit{C}}_{S}{\mathit{C}}_{I}$ and $r}_{PSC}\approx \frac{{\mathit{C}}_{S}{\mathit{C}}_{I}}{{\mathit{C}}_{S}+{\mathit{C}}_{I}$ show the simple dependence of the variance ratio and power spectral centroid ratio on the coupling parameters (visualized in Figure 3I–J). The experimental results, where the mean component had higher variance and lower power spectral centroid than the difference component, correspond to the parameter regime of $0<{\mathit{C}}_{I}<{\mathit{C}}_{S}$. In this regime, consistent with the earlier analysis of the full model, $\mathit{C}}_{I$ acts as positive feedback to the mean component to amplify it and slow it down, and acts negative feedback to the difference component to suppress it and speed it up.
Modeling: Kuramoto model
Request a detailed protocolWe explored an alternative mechanism for interbrain coupling using the Kuramoto model (Strogatz, 2000; Acebrón et al., 2005). In this model, the activity of the two brains are treated as two oscillators, whose phases are coupled:
Here ${\theta}_{1}\left(t\right)$ and ${\theta}_{2}\left(t\right)$ are the phases of the two oscillators at time $t$, and the corresponding activity of the two brains at time $t$ are $\mathit{a}}_{1}\left(t\right)=\frac{sin\left({\theta}_{1}\left(t\right)\right)+1}{2$ and $\mathit{a}}_{2}\left(t\right)=\frac{sin\left({\theta}_{2}\left(t\right)\right)+1}{2$; ${\omega}_{1}$ and ${\omega}_{2}$ are the natural frequencies of the two oscillators; $K$ is the coupling strength between the oscillators: $K>0$ for simulations of onechamber sessions, and $K=0$ for simulations of twochambers sessions.
For each simulation, ${\omega}_{1}$ and ${\omega}_{2}$ were drawn from lognormal distributions with means $\overline{\omega}}_{1$ and $\overline{\omega}}_{2$, respectively, and the same standard deviation ${\sigma}_{\omega}$. $\overline{\omega}}_{1$, $\overline{\omega}}_{2$, and ${\sigma}_{\omega}$ were set to make ${\omega}_{1}$ and ${\omega}_{2}$ different from each other, to avoid the two brains trivially synchronizing by having the same natural frequency. To simulate the Kuramoto model, equation (18) was numerically integrated (ode15s function in MATLAB) with initial condition ${\theta}_{1}\left(0\right)={\theta}_{2}\left(0\right)=0$. The phases ${\theta}_{1}\left(t\right)$ and ${\theta}_{2}\left(t\right)$ were then converted to activity ${\mathit{a}}_{1}\left(t\right)$ and ${\mathit{a}}_{2}\left(t\right)$, which were plotted in Figure 4—figure supplement 1AD and analyzed in Figure 4—figure supplement 1EF. This shows that the phasecoupling mechanism of the Kuramoto model is able to reproduce interbrain correlation and the relative magnitudes of the mean and difference components from the data; however, it does not reproduce the relative timescales of the mean and difference components from the data.
Modeling: correlation between activity variables under rotated activity bases
Request a detailed protocolThe neural dynamics in the model are governed by (2). To express (2) under a different basis, $\mathit{a}\left(t\right)$, $\mathit{C}$, and $\mathit{b}\left(t\right)$ are transformed to ${\mathit{a}}^{\prime}\left(t\right)={\mathit{U}}^{1}\mathit{a}\left(t\right)$, $\mathit{C}}^{\prime}={\mathit{U}}^{1}\mathit{C}\mathit{U$, and ${\mathit{b}}^{\prime}\left(t\right)={\mathit{U}}^{1}\mathit{b}\left(t\right)$, respectively, where the columns of the matrix $\mathit{U}$ are the new basis vectors. If $\mathit{U}=\left(\begin{array}{cc}\text{cos}\theta & \text{sin}\theta \\ \text{sin}\theta & \text{cos}\theta \end{array}\right)$, then the new basis corresponds to a counterclockwise rotation of the original basis (where each axis is the activity of one bat) by an angle of $\theta $. By rotating the basis, the functional coupling between the activity variables (i.e. between the two elements of ${\mathit{a}}^{\prime}\left(t\right)$) changes as a function of the rotation angle (Figure 5A and D). This functional coupling in turn determines the correlations between those activity variables, which forms a set of predictions of the model (Figure 5B and E). These predictions are confirmed in the data in Figure 5C and F.
Additionally, we examined the influence of behavioral modulation on correlation between the activity variables under rotated bases, by regressing out behavior from the activity variables. To do this, for each behavior (e.g. resting, selfgrooming, probing, etc.), each bat, and each simulation or experimental session, we used a binary vector to represent the time course of the bat engaging in that behavior: the elements of the vector correspond to time points, and the values of the elements are ‘1’ at time points when the bat was engaging in that behavior, and ‘0’ when it was not. Then, we performed linear regression to predict each activity variable over time under each basis using the set of all behavioral binary vectors of both bats as predictors. The linear regression fit was then subtracted from the activity, and correlation was calculated between these residuals. The results are shown in Figure 5B, C, E and F (brown lines): for both model and data, regressing out behaviors changes the magnitudes, but not the shapes of the correlation curves.
Modeling: relationship between variance ratio and power spectral centroid ratio
Request a detailed protocolHere we analyze the relationship between the mean/difference ratio of the variance (${r}_{Var}$) and the mean/difference ratio of the power spectral centroid (${r}_{PSC}$) in the model (seen in Figure 5G). ${r}_{Var}$ and ${r}_{PSC}$ are functions of the behavioral modulation $\mathit{b}\left(t\right)$, which differs from simulation to simulation. Figure 5G shows that ${r}_{Var}$ and ${r}_{PSC}$ corresponding to different instances of $\mathit{b}\left(t\right)$ tend to covary linearly. To examine whether this linear relationship is due to a systematic relationship between different instances of $\mathit{b}\left(t\right)$, or due to the neural dynamics, we examine the effect of random perturbations to behavioral modulation on the relationship between ${r}_{Var}$ and ${r}_{PSC}$. To do so, we first express ${r}_{Var}$ and ${r}_{PSC}$ in the frequency domain, and then examine whether and how they covary given random perturbations.
We first follow the procedure of the earlier section on the reduced model to approximate the variance and the power spectral centroid. In that section, we analyzed the reduced model with $\tau =1$ and ${\left{\stackrel{\sim}{\mathit{b}}}_{M,k}\right}^{2}={\left{\stackrel{\sim}{\mathit{b}}}_{D,k}\right}^{2}=p$ for $k\ge 1$; here we consider arbitrary $\tau$, $\left{\stackrel{\sim}{\mathit{b}}}_{M,k}\right}^{2$, and $\left{\stackrel{\sim}{\mathit{b}}}_{D,k}\right}^{2$. In this case, the neural decay functions are $f}_{M}\left(t\right)={e}^{\frac{{\mathit{C}}_{S}{\mathit{C}}_{I}}{\tau}t$ and $f}_{D}\left(t\right)={e}^{\frac{{\mathit{C}}_{S}+{\mathit{C}}_{I}}{\tau}t$, and the variance and power spectral centroid of the mean component are
where $\mathit{B}}_{M,k}={\left{\stackrel{\sim}{\mathit{b}}}_{M,k}\right}^{2$ and ${F}_{M,k}={\left{\stackrel{~}{f}}_{M,k}\right}^{2}$. The expressions for the difference component are similar.
Then, the mean/difference ratio of the variance is
The mean/difference ratio of the power spectral centroid is
To examine how ${r}_{Var}$ and ${r}_{PSC}$ change with changes in behavioral modulation, we calculate the following partial derivatives. The partial derivative of ${r}_{Var}$ with respect to the power of the mean component of behavioral modulation at frequency $k$ is
The partial derivative of ${r}_{Var}$ with respect to the power of the difference component of behavioral modulation at frequency $k$ is
The partial derivative of ${r}_{PSC}$ with respect to the power of the mean component of behavioral modulation at frequency $k$ is
The partial derivative of ${r}_{PSC}$ with respect to the power of the difference component of behavioral modulation at frequency $k$ is
In subsequent calculations, the infinite sums were truncated at $k}_{t$ such that $\frac{{k}_{t}}{T}=0.2$ Hz, which is the Nyquist frequency for the sampled simulated activity used for our analyses.
We now consider perturbations to the power spectra of behavioral modulation. We concatenate the power spectra of the behavioral modulation of the mean and difference components to form
We perturb a given $\mathit{B}$ by adding $\delta \mathit{B}$, drawn from a uniform distribution on the hypersphere centered at 0. To see whether and how ${r}_{Var}$ and ${r}_{PSC}$ covary given random perturbations to behavioral modulation, we next determine the covariance matrix between the changes in ${r}_{Var}$ and ${r}_{PSC}$ resulting from perturbations $\delta \mathit{B}$.
Considering ${r}_{Var}$ and ${r}_{PSC}$ as functions of $\mathit{B}$, their gradients are
where the partial derivatives are given in (21)(24). Given a small perturbation $\delta \mathit{B}$, the resulting changes to ${r}_{Var}$ and ${r}_{PSC}$ are approximately $\delta {\mathit{B}}^{T}\mathrm{\nabla}{r}_{Var}$ and $\delta {\mathit{B}}^{T}\mathrm{\nabla}{r}_{PSC}$, respectively. To calculate the covariance matrix between $\delta {\mathit{B}}^{T}\mathrm{\nabla}{r}_{Var}$ and $\delta {\mathit{B}}^{T}\mathrm{\nabla}{r}_{PSC}$, we first perform an orthogonal transformation to any orthonormal basis whose first basis vector is $\nabla {r}_{Var}/\Vert \nabla {r}_{Var}\Vert $, where $\Vert \bullet \Vert $ denotes vector norm. We use $\delta \hat{\mathit{B}}$, $\mathrm{\nabla}{\hat{r}}_{Var}$, and $\mathrm{\nabla}{\hat{r}}_{PSC}$ to denote $\delta \mathit{B}$, $\nabla {r}_{Var}$, and $\nabla {r}_{PSC}$ under the new basis, respectively. The variance of $\delta {\mathit{B}}^{T}\mathrm{\nabla}{r}_{Var}$ is
Similarly, the variance of $\delta {\mathit{B}}^{T}\mathrm{\nabla}{r}_{PSC}$ is $\frac{{\Vert \delta \mathit{B}\Vert}^{2}}{2{k}_{t}}{\Vert \mathrm{\nabla}{r}_{PSC}\Vert}^{2}$. The covariance between $\delta {\mathit{B}}^{T}\mathrm{\nabla}{r}_{Var}$ and $\delta {\mathit{B}}^{T}\mathrm{\nabla}{r}_{PSC}$ is
Thus, the covariance matrix between $\delta {\mathit{B}}^{T}\mathrm{\nabla}{r}_{Var}$ and $\delta {\mathit{B}}^{T}\mathrm{\nabla}{r}_{PSC}$ is
The eigenvector of this matrix corresponding to the larger eigenvalue is the direction of the local linear trend, whereas the relative sizes of the eigenvalues indicate the strength of the linear trend (the larger the difference between the eigenvalues, the stronger the linear trend).
We examined the eigenvectors and eigenvalues of (25) across different simulations (Figure 5—figure supplement 1; in computing (25), integrals were evaluated numerically). We found that ${r}_{Var}$ and ${r}_{PSC}$ resulting from random perturbations to $\mathit{B}$ consistently show linear relationships. Furthermore, the slopes of these local linear relationships, which are not influenced by systematic variations of behavioral modulation $\mathit{b}\left(t\right)$ across simulations, are consistent with the slopes of the global linear relationships seen in Figure 5G.
Modeling: using the neural activity of one bat to discriminate the behavior of the other bat
Request a detailed protocolIn the model, the neural activity of each bat is directly modulated by its own behavior (e.g. ${\mathit{a}}_{1}\left(t\right)$ is modulated by ${\mathit{b}}_{1}\left(t\right)$). Additionally, in the onechamber model, the activity of each bat is indirectly modulated by the behavior of the other bat, through functional acrossbrain coupling (e.g. ${\mathit{a}}_{1}\left(t\right)$ is modulated by ${\mathit{b}}_{2}\left(t\right)$ through the coupling $\mathit{C}}_{I$). This naturally suggests that the neural activity of each bat should encode the behavior of the other bat independently from encoding its own behavior.
We used the following method to quantify this. Before examining whether the activity of one bat can discriminate the behavior of the other bat, we first regressed out the behavior of each bat from its own neural activity, using the method described in the earlier section ‘Modeling: correlation between activity variables under rotated activity bases’. Importantly, this eliminates potential spurious correlation between one bat’s activity and another bat’s behavior caused by any coordinated behaviors between the bats. We then asked whether the neural activity of a given bat discriminates between a given pair of behaviors by the other bat (e.g., using bat 1’s neural activity to discriminate whether bat 2 is resting or engaging in social grooming), by plotting the receiver operating characteristic (ROC) curve and calculating the area under the curve (Figure 6B and E; Dayan and Abbott, 2005). For each discrimination, positive and negative class assignments for the two behaviors were made so that the area under the ROC curve is greater than or equal to 0.5. We then averaged the area under the ROC curve across bats, pairs of behaviors, and simulations or sessions (Figure 6C and F). This showed that the activity of each bat was modulated by the behaviors of the other bat independently of its own behavior, in both the model and the data.
We note that for both the model and the data, the area under the ROC curve was above 0.5 for the twochambers sessions—this is a necessary consequence of noise and finite sample size. To illustrate this, we can consider a hypothetical example where neural activity does not encode behaviors A and B, that is, the distributions of neural activity during these two behaviors are identical. If we have infinite amount of data, or if the neural activity is noiseless during these two behaviors, then the empirically observed distributions of activity during these two behaviors would be exactly identical, so that the area under the ROC curve would be exactly 0.5. However, in reality we have finite amounts of noisy data, so that the two empirically observed distributions would not be identical. These different empirical distributions would then necessarily result in an area under the ROC curve greater than 0.5, since positive and negative classes were assigned to the two behaviors such that the area under the ROC curve is greater than or equal to 0.5.
Modeling: interbrain relationship during group social interactions
Request a detailed protocolTo generalize our twobat model to more than two bats, we used the same equation (2), with $\mathit{a}\left(t\right)$ and $\mathit{b}\left(t\right)$ now being $n$dimensional vectors, and $\mathit{C}$ being an $n\times n$ matrix, where $n$ is the number of interacting bats. $\mathit{C}$ retains the same structure from the twobat model: all diagonal elements (functional selfcoupling) are ${\mathit{C}}_{S}$, and all offdiagonal elements (functional acrossbrain coupling) are $\mathit{C}}_{I$.
To understand the $n$bat model, we examine the eigenvectors and eigenvalues of $\mathit{C}$. Note that the $n$dimensional vector whose elements are all 1s is an eigenvector, with eigenvalue $\left(n1\right){\mathit{C}}_{I}{\mathit{C}}_{S}$. This eigenvector corresponds to the direction of the mean activity across all bats, and is thus the $n$bat analogue of the mean component from the twobat model. Any vector orthogonal to the $n$bat mean component is also an eigenvector, with eigenvalue ${\mathit{C}}_{I}{\mathit{C}}_{S}$. These eigenvectors define an $\left(n1\right)$dimensional subspace, which contains all interbrain activity patterns that correspond to activity differences across brains; we call this subspace the difference subspace, which is the multidimensional analogue of the difference component from the twobat model. To ensure stability (so that neural activity do not go to infinity), we take our parameter regime to be $0<\left(n1\right){\mathit{C}}_{I}<{\mathit{C}}_{S}$. Thus, ${\mathit{C}}_{I}{\mathit{C}}_{S}<\left(n1\right){\mathit{C}}_{I}{\mathit{C}}_{S}<0$. This means that, similar to the twobat model, the $n$bat mean component is amplified and slowed down by the positive feedback provided by functional acrossbrain coupling, whereas activity patterns in the difference subspace are suppressed and sped up by negative feedback. This results in the predictions that, for $n$ socially interacting bats, activity in the $n$bat mean component will have larger magnitude and slower timescales than activity in the difference subspace on average (Figure 7C–D).
Another prediction concerns the correlation between activity variables during group interactions. Because the activity of pairs of bats are positively functionally coupled in the $n$bat model, the model would predict positive interbrain correlations (Figure 7E). On the other hand, because the $n$bat mean component and vectors in the difference subspace are all eigenvectors, they are functionally uncoupled. Thus, the model would predict zero average correlation between the $n$bat mean component and activity variables in the difference subspace (Figure 7E).
For the simulations in Figure 7C–E, we used an $n$bat model with $n=4$ and performed the simulations using the same procedures as for the twobat model. Because we do not have behavioral statistics for fourbat group interactions, we opted to use noise for $\mathit{b}\left(t\right)$, so that model predictions based on functional acrossbrain coupling would not be biased by structures in $\mathit{b}\left(t\right)$ (e.g. the $n$bat mean component could have larger magnitudes than the difference subspace if $\mathit{b}\left(t\right)$ has such a structure). To generate $\mathit{b}\left(t\right)$ for the $n$bat model, we used the following procedure. For a given simulated session and a given bat i, we generated ${\mathit{b}}_{i}\left(t\right)$ independently of the other bats. We generated a random vector $\mathit{b}}_{pre$, whose dimensionality was the number of time points spaced 2.5 s apart in the simulated session (matching the time step size of the Markov chain from the twobat model). The elements of $\mathit{b}}_{pre$ were drawn independently from a Gaussian distribution with mean $b}_{mean$ and standard deviation $b}_{std$. We smoothed $\mathit{b}}_{pre$ with a 1200point moving average filter, then added independent Gaussian noise to its elements (0 mean, standard deviation ${\sigma}_{n}$) as in the twobat model. The resulting vector contained the values of ${\mathit{b}}_{i}\left(t\right)$ at 2.5 s intervals; ${\mathit{b}}_{i}\left(t\right)$ at time points inbetween were linearly interpolated. The same procedure was repeated for each bat. Note that, because $\mathit{b}\left(t\right)$ was randomly generated noise, the simulations only offered qualitative predictions. Indeed, the data confirmed the qualitative trends seen in the simulations, but differed from them quantitatively (Figure 7C–H).
Modeling: parameters
Request a detailed protocolIn our data, all four neural signals (30–150 Hz LFP power, 1–29 Hz LFP power, multiunits, and single units) show the same qualitative phenomena: faster and smaller difference component compared to the mean component, on onechamber sessions but not twochambers sessions (Figure 1). The four signals differed quantitatively in the extent they showed the same qualitative phenomena. The goal of our model is not to quantitatively fit any one of the four neural signals in particular, but to provide mechanistic explanations for the general qualitative phenomena that do not require finetuning of parameters. As explained above, the opposite feedback mechanism depends simply on $0<{\mathit{C}}_{I}<{\mathit{C}}_{S}$ , whereas the coordinated behavioral modulation mechanism is a manifestation of the empirical behavioral transition frequencies.
The focus of the model is to reproduce and explain the qualitative trend of the difference component being faster and smaller than the mean component. Other aspects of the data could be reproduced by extending the model with additional parameters, but we chose not to do so in the interest of keeping the model simple and focused. For example, the 30–150 Hz LFP power data showed more variability across sessions compared to the model simulations (Figure 3E–F). The higher variability could be reproduced if we introduced variability across simulations in the behavioral transition probabilities or the strength of functional acrossbrain coupling.
The model parameters for the twobat models are: ${\mathit{C}}_{S}=1$, ${\mathit{C}}_{I}=0.4$ for simulations with functional acrossbrain coupling, ${\mathit{C}}_{I}=0$ for simulations without functional acrossbrain coupling, $\tau =15s$, ${\sigma}_{n}=0.15$, $T=100$ minutes, ${b}_{resting}=0.158$, ${b}_{active\text{}nonsocial}=0.269$, ${b}_{selfgrooming}=0.264$, ${b}_{social\text{}grooming}=0.223$, ${b}_{probing}=0.284$, ${b}_{fighting}=0.355$, ${b}_{mating}=0.367$, ${b}_{wing\text{}covering}=0.364$, ${b}_{reaching}=0.321$, ${b}_{blocking}=0.339$, ${b}_{other\text{}interactions}=0.291$, ${b}_{constant}=0.08$ for onechamber simulations and ${b}_{constant}=0$ for twochambers simulations. The behavioral modulation parameter for each behavior listed above (${b}_{resting}$, etc.) was set as the average 30–150 Hz mean normalized LFP power during that behavior from the data: take a given bat, for each session, average its 30–150 Hz mean normalized LFP power across all channels, pool together these averaged power values from all time points when the bat engaged in the given behavior from all sessions (including both onechamber and twochambers), then pool across all bats, and then average across all such pooled data. ${b}_{constant}$, $\mathit{C}}_{S$, $\mathit{C}}_{I$, $\tau $, and ${\sigma}_{n}$ were chosen so that the levels of simulated neural activity roughly match the levels of 30–150 Hz mean normalized LFP power observed during the experiments.
The model parameters for the fourbat model are: ${\mathit{C}}_{S}=1$, ${\mathit{C}}_{I}=0.1$, ${b}_{mean}=0.2$, and ${b}_{std}=3.5$. All other parameters are the same as in the twobat model.
The model parameters for the Kuramoto model are: $K=0.0035$ for simulations of onechamber sessions, $K=0$ for simulations of twochambers sessions, ${\overline{\omega}}_{1}=0.005$, ${\overline{\omega}}_{2}=0.01$, and ${\sigma}_{\omega}=0.0005$.
Sample sizes
Request a detailed protocolIn this section we list the sample sizes for all results that were presented as averages.
Figure 1M, N, Q and R: n = 52 sessions for onechamber sessions, and 18 sessions for twochambers sessions.
Figure 1O and S: n = 675 multiunit pairs for onechamber sessions, and 284 multiunit pairs for twochambers sessions.
Figure 1P and T: n = 256 single unit pairs for onechamber sessions, and 65 single unit pairs for twochambers sessions.
Figure 3E–F: for data, n = 52 sessions for onechamber sessions, and 18 sessions for twochambers sessions; for model, n = 100 simulations for onechamber model, and 100 simulations for twochambers model.
Figure 4G: for data, n = 50 sessions; for model, n = 100 simulations with functional acrossbrain coupling, and 100 simulations without functional acrossbrain coupling.
Figure 5B: n = 100 simulations.
Figure 5C: n = 52 sessions.
Figure 5E: n = 100 simulations.
Figure 5F: n = 18 sessions.
Figure 6C: n = 8,235 simulations × bats × behavior pairs for onechamber simulations, n = 1,848 simulations × bats × behavior pairs for twochambers simulations.
Figure 6F: n = 2086 sessions × bats × behavior pairs for onechamber sessions, n = 200 sessions × bats × behavior pairs for twochambers sessions.
Figure 7C–E: n = 100 simulations.
Figure 7F–H: n = 20 sessions.
Figure 3—figure supplement 2A, B: n = 100 crossvalidation test sets.
Figure 3—figure supplement 2H, I: for data, n = 52 sessions for onechamber sessions, and 18 sessions for twochambers sessions; for model, n = 100 simulations for onechamber model, and 100 simulations for twochambers model.
Figure 4—figure supplement 1E, F: for data, n = 52 sessions for onechamber sessions, and 18 sessions for twochambers sessions; for the Kuramoto model, n = 100 simulations for onechamber model, and 100 simulations for twochambers model.
Figure 5—figure supplement 1C, E: n = 100 simulations for onechamber model, and 100 simulations for twochambers model.
Data availability
Source code for the models is available at https://github.com/zhangwujie/Neurobatlabcodes/tree/master/Interbrainmodel, copy archived at swh:1:rev:20f445028b5b1fa3b8b15be020252a397eb6479f.
References

The Kuramoto model: A simple paradigm for synchronization phenomenaReviews of Modern Physics 77:137–185.https://doi.org/10.1103/RevModPhys.77.137

The neurobiology of social cognitionCurrent Opinion in Neurobiology 11:231–239.https://doi.org/10.1016/s09594388(00)002026

Meeting of minds: the medial frontal cortex and social cognitionNature Reviews. Neuroscience 7:268–277.https://doi.org/10.1038/nrn1884

Social neuroscience and hyperscanning techniques: past, present and futureNeuroscience and Biobehavioral Reviews 44:76–93.https://doi.org/10.1016/j.neubiorev.2012.07.006

Generative models of cortical oscillations: neurobiological implications of the kuramoto modelFrontiers in Human Neuroscience 4:1–14.https://doi.org/10.3389/fnhum.2010.00190

Neuronal reference frames for social decisions in primate frontal cortexNature Neuroscience 16:243–250.https://doi.org/10.1038/nn.3287

A taxonomy of external and internal attentionAnnual Review of Psychology 62:73–101.https://doi.org/10.1146/annurev.psych.093008.100427

A selective review of selective attention research from the past centuryBritish Journal of Psychology (London, England 92:53–78.https://doi.org/10.1348/000712601162103

Resource Ephemerality Drives Social Foraging in BatsCurrent Biology 28:3667–3673.https://doi.org/10.1016/j.cub.2018.09.064

Marmoset vocal communication: Behavior and neurobiologyDevelopmental Neurobiology 77:286–299.https://doi.org/10.1002/dneu.22464

The role of the human prefrontal cortex in social cognition and moral judgmentAnnual Review of Neuroscience 33:299–324.https://doi.org/10.1146/annurevneuro060909153230

Social interaction networks in the primate brainCurrent Opinion in Neurobiology 65:49–58.https://doi.org/10.1016/j.conb.2020.08.012

Auditory attentionfocusing the searchlight on soundCurrent Opinion in Neurobiology 17:437–455.https://doi.org/10.1016/j.conb.2007.07.011

Global forebrain dynamics predict rat behavioral states and their transitionsThe Journal of Neuroscience 24:11137–11147.https://doi.org/10.1523/JNEUROSCI.352404.2004

Persistent producerscrounger relationships in batsScience Advances 4:e1603293.https://doi.org/10.1126/sciadv.1603293

Braintobrain coupling: a mechanism for creating and sharing a social worldTrends in Cognitive Sciences 16:114–121.https://doi.org/10.1016/j.tics.2011.12.007

Mirroring and beyond: coupled dynamics as a generalized framework for modelling social interactionsPhilosophical Transactions of the Royal Society of London. Series B, Biological Sciences 371:20150366.https://doi.org/10.1098/rstb.2015.0366

Getting to know you: reputation and trust in a twoperson economic exchangeScience (New York, N.Y.) 308:78–83.https://doi.org/10.1126/science.1108062

A MultiBrain Framework for Social InteractionTrends in Neurosciences 43:651–666.https://doi.org/10.1016/j.tins.2020.06.008

The twobrain approach: how can mutually interacting brains teach us something about social interaction?Frontiers in Human Neuroscience 6:1–10.https://doi.org/10.3389/fnhum.2012.00215

Perception of social synchrony induces motherchild gamma coupling in the social brainSocial Cognitive and Affective Neuroscience 12:1036–1046.https://doi.org/10.1093/scan/nsx032

DeepLabCut: markerless pose estimation of userdefined body parts with deep learningNature Neuroscience 21:1281–1289.https://doi.org/10.1038/s415930180209y

Responses of primate frontal cortex neurons during natural vocal communicationJournal of Neurophysiology 114:1158–1171.https://doi.org/10.1152/jn.01003.2014

The role of gamma interbrain synchrony in social coordination when humans face territorial threatsSocial Cognitive and Affective Neuroscience 12:1614–1623.https://doi.org/10.1093/scan/nsx093

Social ContextDependent Activity in Marmoset Frontal Cortex Populations during Natural ConversationsThe Journal of Neuroscience 37:7036–7047.https://doi.org/10.1523/JNEUROSCI.070217.2017

Social placecells in the bat hippocampusScience (New York, N.Y.) 359:218–224.https://doi.org/10.1126/science.aao3474

Neuronal correlates of strategic cooperation in monkeysNature Neuroscience 24:116–128.https://doi.org/10.1038/s41593020007469

Fast animal pose estimation using deep neural networksNature Methods 16:117–125.https://doi.org/10.1038/s4159201802345

Infant and Adult Brains Are Coupled to the Dynamics of Natural CommunicationPsychological Science 31:6–17.https://doi.org/10.1177/0956797619878698

Using secondperson neuroscience to elucidate the mechanisms of social interactionNature Reviews. Neuroscience 20:495–505.https://doi.org/10.1038/s4158301901794

Attentional modulation of visual processingAnnual Review of Neuroscience 27:611–647.https://doi.org/10.1146/annurev.neuro.26.041002.131039

Cortical representation of group social communication in batsScience (New York, N.Y.) 374:eaba9584.https://doi.org/10.1126/science.aba9584

The contribution of distinct subregions of the ventromedial frontal cortex to emotion, social behavior, and decision makingCognitive, Affective & Behavioral Neuroscience 8:485–497.https://doi.org/10.3758/CABN.8.4.485

What can we learn from a twobrain approach to verbal interaction?Neuroscience and Biobehavioral Reviews 68:454–459.https://doi.org/10.1016/j.neubiorev.2016.06.009

Interindividual synchronization of brain activity during live verbal communicationBehavioural Brain Research 258:75–79.https://doi.org/10.1016/j.bbr.2013.10.015

From the field to the lab and back: neuroethology of primate social behaviorCurrent Opinion in Neurobiology 68:76–83.https://doi.org/10.1016/j.conb.2021.01.005

Agentspecific responses in the cingulate cortex during economic exchangesScience (New York, N.Y.) 312:1047–1050.https://doi.org/10.1126/science.1125596

Social DecisionMaking and the Brain: A Comparative PerspectiveTrends in Cognitive Sciences 21:265–276.https://doi.org/10.1016/j.tics.2017.01.007

Interpersonal Neural Entrainment during Early Social InteractionTrends in Cognitive Sciences 24:329–342.https://doi.org/10.1016/j.tics.2020.01.006

Representation of threedimensional space in the hippocampus of flying batsScience (New York, N.Y.) 340:367–372.https://doi.org/10.1126/science.1235338

How we transmit memories to other brains: constructing shared neural representations via communicationCerebral Cortex (New York, N.Y.: 1991) 27:4988–5000.https://doi.org/10.1093/cercor/bhx202

History of winning remodels thalamoPFC circuit to reinforce social dominanceScience (New York, N.Y.) 357:162–168.https://doi.org/10.1126/science.aak9726
Decision letter

Noah J CowanReviewing Editor; Johns Hopkins University, United States

Ronald L CalabreseSenior Editor; Emory University, United States
In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.
Decision letter after peer review:
Thank you for submitting your article "A unifying mechanism governing interbrain neural relationship during social interactions" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, including Noah Cowan as the Reviewing Editor and Reviewer #1, and the evaluation has been overseen by Ronald Calabrese as the Senior Editor. Two of the reviewers have opted to remain anonymous.
The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this decision letter.
After careful consideration and discussion, the reviewers concur that there is some potential for this work, but the reviewers are unanimous in noting that the paper has serious weaknesses and requires extensive new analysis and a major revision to be suitable for publication in eLife. We note that this in part stems from the high bar of publication in eLife; this paper may be suitable in a more field specific journal as is, and so the authors will need to balance the requirements we are placing for new revisions and the desire to publish the work as is.
Also, it was noted that this manuscript reuses portions of the methods section of a prior article by the same authors (Zhang & Yartsev 2019 Cell). The methods should be rewritten so as to avoid this and it should be made clear when you are quoting directly from the prior manuscript.
If the authors choose to revise the manuscript for eLife, it is not guaranteed that future revisions will be accepted, as the model would still need to be vetted by the reviewers.
Reviewer #1:
This paper provides experimental and modeling analysis of the interbrain coupling of socially interacting bats, and reports that coordinated brain activity evolves at a slower time scale than the activity describing the differences. Specifically, the paper finds that there is an attracting submanifold corresponding to the mean (or "common mode") of neural activity, and that the dynamics in the orthogonal eigenmode, corresponding to the difference in brain activity, decays rapidly. These rapid decays in the difference mode are referred to as "catch up" activity.
There are two main findings:
1) Neural activity (especially higher frequency LFP activity in the 30150Hz range) is modulated by social context. Specifically, the ratio of the averaged, momenttomoment MEAN:DIFF ratio is much higher when the bats are in a single chamber, clearly indicating that the animals are coordinating their neural activity. This change also seems to hold  although not as striking  in lowerfrequency LFP and spiking activity.
2) The time scales of the mean vs. difference dynamics are segregated: the "difference dynamics" evolve at a faster time scale than "similarity dynamics", seems to be well supported.
The basic finding is presented in Figure 1. The rest of the paper is focused on a modeling study to garner further insight into the dynamics.
Weaknesses:
This is an entirely phenomenological paper, and while it claims to garner "mechanistic insight", it is unclear what that means.
The basic idea of the model is simple and somewhat interesting, but the details are extremely complex. There are many examples of this, but the method used to "regress out" the behavior was very hard to interpret.
On the face of it, the model is extremely simple: a twostate linear dynamical system. However, this simplistic description buries extreme complexity. The model is extremely complex as involves a large number of parameters (e.g., time switching 'b' values, the values of which are completely unclear), the switching over time of these parameters based on handscored animal behavioral state, and the complex mix of markovian and linear dynamical systems theoretic results. Indeed, a fundamental weakness of the model is that the Markov chain is taken as an "input" to the 2state linear systems model, as if somehow the neural state does not affect the state transitions. Further, the Markov assumption is not rigorously tested. No model selecting or other model validation appears to be done.
In short, the model, while very interesting, is so complex that it is literally impossible to evaluate. The authors report literally no shortcomings of their model. They do not report parameter estimation methods. They do not report fitting errors or other model validation metrics. The only evaluation is whether it can produce certain outputs that are similar to biological data. While the latter is certainly important, all models are wrong, and it essential to have a model simple enough to understand, both in terms of how it works and how it fails.
In general, while the basic finding is fairly interesting, and the experiments and their findings are highly relevant to the field, the modeling and its explication fall short.
It is not that it is wrong or bad; however, it is not clear that such a complex model increases our understanding beyond the experimental findings in Figure 1, and if it does, there has to be a major caveat that the model itself is not carefully vetted.
Reviewer #2:
In this paper, Wujie Zhang and Michael Yartsev investigate some of the basic underpinnings of interbrain synchrony in socially interacting animals. The phenomena of inter brain synchrony is fascinating and has been observed in a variety of situations across different mammalian species. It has also been proposed to play a critical role in certain social behaviors. Here, the authors report that brain activity across two interacting bats display not only similarities but also important differences. The also use advance computational modeling to capture s these two characteristics as well as to demonstrate how they are affected by the presence and absence of interaction between animal pairs.
Reviewer #3:
The activity in the frontal cortex of mammals has been previously shown to become more correlated in socially interacting animals than when they are alone. In the current study, the authors examine the differences in brain activity that emerge during social interactions. The correlations and differences in activity were shown to occur over different time scales, with mean correlations occurring over longer time scales whereas differences occur over shorter time scales. The authors made a model of these processes that show how feedback may give rise to these phenomena.
https://doi.org/10.7554/eLife.70493.sa1Author response
Essential revisions:
1) Extensive rewrite of the manuscript to improve clarity.
2) Extensive analysis of the modeling.
3) Clearer distinction from the prior work (including not copying prior methods) by the same group.
We have revised our manuscript extensively, based on the three points of essential revisions above and following the reviewers’ suggestions, as detailed in our pointbypoint reply. Furthermore, we have added a new data set of neural recording from a larger social group (four bats) engaged in free social interactions, which confirmed the predictions of our model.
Please see the individual reviewer comments for guidance on how to rewrite the manuscript.
As for Item #2, Reviewer #1 raised most of the concerns, and they were discussed among all of the reviewers who concur. During the discussion, it was raised that perhaps the findings of the model didn't really increase understanding past the basic data finding in Figure 1. In particular, there are a huge number of fitting parameters, such as the "b" values. The b values seem to effectively change the DC set point of the model for each handscored behavioral epoch, via
a* =  C^{}}_{b}
(this is not explained in the manuscript). This is roughly equivalent to taking all events ("grooming", "mating", etc), finding nominal DC point of the associated LFP, and subtracting it off, i.e. the b terms add an affine term that biases the linear dynamics to a nominal value. This makes sense, since power is a nonnegative quantity, and adding this DC offset effectively linearizes an affine system around a new operating point for each behavioral epoch, allowing for linear regression on the time constants.
We could not determine how the authors fit the model parameters. If the parameters were optimized to reproduce the predictions, as opposed to fit the time domain data, for example, the model would be unconvincing and unfortunately not suitable for eLife; there are simply too many model DOFs and too few terms to predict, and so it is a foregone conclusion that you can fit the outcome.
This speaks to the broader point that there was agreement among the reviewers that there was no model validation. For example, we could not find reported even a single timedomain trace comparing the model with data.
Again, how does the model move past Figure 1 data analysis? Instead of fitting a complex, hierarchical model, would it be sufficient to analyze the constants by looking at the autocorrelation functions for each term (a2 + a1)/2 and (a2a1)/2? The paper basically does this in terms of the PSD (Fourier transform of the autocorrelation functions).
Reviewer #1:
We thank the reviewer for the helpful comments and constructive criticism, which has helped us better understand our model and improve our manuscript.
This paper provides experimental and modeling analysis of the interbrain coupling of socially interacting bats, and reports that coordinated brain activity evolves at a slower time scale than the activity describing the differences. Specifically, the paper finds that there is an attracting submanifold corresponding to the mean (or "common mode") of neural activity, and that the dynamics in the orthogonal eigenmode, corresponding to the difference in brain activity, decays rapidly. These rapid decays in the difference mode are referred to as "catch up" activity.
There are two main findings:
1) Neural activity (especially higher frequency LFP activity in the 30150Hz range) is modulated by social context. Specifically, the ratio of the averaged, momenttomoment MEAN:DIFF ratio is much higher when the bats are in a single chamber, clearly indicating that the animals are coordinating their neural activity. This change also seems to hold  although not as striking  in lowerfrequency LFP and spiking activity.
2) The time scales of the mean vs. difference dynamics are segregated: the "difference dynamics" evolve at a faster time scale than "similarity dynamics", seems to be well supported.
The basic finding is presented in Figure 1. The rest of the paper is focused on a modeling study to garner further insight into the dynamics.
Weaknesses:
This is an entirely phenomenological paper, and while it claims to garner "mechanistic insight", it is unclear what that means.
We regret not clarifying sufficiently what we meant by “mechanistic insight.” The insight is the following: functional acrossbrain coupling acts as positive feedback to the mean component of neural activity, which amplifies it and slows it down; at the same time, it acts as negative feedback to the difference component, which suppresses it and speeds it up. Thus, findings (1) and (2) in the reviewer’s summary above can be explained by the same model mechanism.
As the reviewer pointed out below, the details of the model are complex, which could have made the simple mechanism above opaque. Thus, we analyzed two simplified versions of the model to make the mechanistic insight clear. This is detailed below in our response to the reviewer’s comment on model complexity.
The basic idea of the model is simple and somewhat interesting, but the details are extremely complex. There are many examples of this, but the method used to "regress out" the behavior was very hard to interpret.
The method for regressing out behavior was described in Materials and methods, and we regret having neglected to reference it in the main text. We now reference it at the first instance in the main text where this is relevant.
On the face of it, the model is extremely simple: a twostate linear dynamical system. However, this simplistic description buries extreme complexity. The model is extremely complex as involves a large number of parameters (e.g., time switching 'b' values, the values of which are completely unclear), the switching over time of these parameters based on handscored animal behavioral state, and the complex mix of markovian and linear dynamical systems theoretic results.
As the reviewer pointed out, the core of the model is very simple: a linear dynamical system that models neural activity coupling. The model mechanism of positive and negative feedback, which is responsible for reproducing the two experimental results summarized by the reviewer above, is contained in this core (see Materials and methods for details). On top of this, the model has a layer of complexity, involving a Markov chain model of behavior and a large number of behavioral parameters. This layer of complexity is independent from the feedback mechanism of the core of the model. Thus, while it makes the model more biologically realistic, it is not required to reproduce the two main experimental results. To explicitly show this, and to better understand the dependence of model behavior on its parameters, we analyzed two reduced versions of the model.
The first reduced model replaces the behavioral inputs with white noise. The original model is $\tau \frac{d\mathit{a}}{dt}=\mathit{C}\mathit{a}+\mathit{b}$, where $\mathit{a}$ is neural activity, $\mathit{C}=\left(\begin{array}{cc}{\mathit{C}}_{S}& {\mathit{C}}_{I}\\ {\mathit{C}}_{I}& {\mathit{C}}_{S}\end{array}\right)$ is the coupling matrix, $\mathit{b}$ is behavioral modulation, and $\tau$ is a time constant. $\mathit{b}$ is where the complexity lies, as it is simulated using a Markov chain and involves many parameters. To strip away this layer of complexity, we replaced $\mathit{b}$ with noise having a simple structure, namely, the mean and difference components of $\mathit{b}$ having identical, flat power spectra. Importantly, this noise input does not induce correlation between bats, and it amounts to inputs of the same magnitude and same timescales to the mean and difference components of $\mathit{a}$.
The resulting reduced model has only two parameters, the functional selfcoupling $\mathit{C}}_{S$ and functional acrossbrain coupling $\mathit{C}}_{I$ (for simplicity, $\tau$ can be absorbed into the other parameters). We are interested in the two results the reviewer summarized above: (1) the mean component of neural activity having a larger variance than the difference component; (2) the mean component having a slower timescale than the difference component. In the manuscript, these are respectively quantified using the variance ratio and the power spectral centroid ratio of the mean and difference components. The reduced model allowed us to derive analytical expressions for these two quantities (see Materials and methods for details). We found that they have very simple dependence on the functional coupling parameters: the variance ratio (mean variance divided by difference variance) is approximately $\frac{{\mathit{C}}_{S}+{\mathit{C}}_{I}}{{\mathit{C}}_{S}{\mathit{C}}_{I}}$, and the centroid ratio (mean centroid divided by difference centroid) is approximately $\frac{{\mathit{C}}_{S}{\mathit{C}}_{I}}{{\mathit{C}}_{S}+{\mathit{C}}_{I}}$.
This parameter dependence is visualized in Figure 3IJ (note that the color maps are in log scale, and the white spaces are regions where the model is unstable). In the experimental data, the mean component had larger variance and lower power spectral centroid than the difference component. This corresponds to the parameter regime of $0<{\mathit{C}}_{I}<{\mathit{C}}_{S}$ (enclosed by dashed lines in Figure 3IJ). Thus, a positive $\mathit{C}}_{I$ acts as positive feedback to the mean component and negative feedback to the difference component, modulating their variance and timescales in opposite directions. This is consistent with the analysis of the original model in Materials and methods. In the revised manuscript, we’ve now added analysis of this reduced model to the Results section.
The reviewer has stated a concern regarding the large number of parameters that set the input level according to behavioral state ($b}_{resting$, $b}_{social\text{}grooming$, $b}_{fighting$, etc.). These parameters are important for ensuring that the model outputs realistic levels of behaviorally modulated neural activity (discussed below in our reply regarding model fit), but they are not important for the main results on variance and timescales. To demonstrate this, we studied a second reduced model. This model is identical to our original model except that, for each simulation, each of the behavior parameters ($b}_{fighting$, etc.) was independently drawn from the uniform distribution from 0 to 1. Despite the completely random behavioral parameters, this reduced model reproduces the variance and timescales results just like the original model, as shown in Author response image 1 (compare with Figure 3EF).
To summarize, the reduced models allowed us to identify the simple parameter dependence of the modeling results, and showed that the simple linear dynamical system at the core of the original model is sufficient to reproduce the two main experimental observations.
Indeed, a fundamental weakness of the model is that the Markov chain is taken as an "input" to the 2state linear systems model, as if somehow the neural state does not affect the state transitions.
Yes, this is a limitation of our model. We’ve added a discussion of this limitation, as well as future directions for overcoming it, in the Discussion section.
The reason we did not model neural control of behavioral transitions is that it is underconstrained by existing data. While the brain obviously controls behaviors, not every part of the brain controls every behavior. Of the 11 behaviors observed in this study, we do not know which of them is controlled by the bat frontal cortex, and we do not know how they might be controlled (i.e., what specific spatiotemporal activity patterns affects behaviors in what ways). Without this knowledge, it’s unclear how to implement neural control of behavior in the model. This knowledge requires perturbation studies (lesion, inactivation, or activity manipulation) to establish casual relationships from neural activity to specific behaviors in the bat, which will be an important future direction.
On the other hand, as the reviewer stated, our model included behavioral modulation of neural activity. It is well known that in mammals, arousal and movement modulate neural activity globally across cortex (McGinley et al., 2015, Neuron). Thus, given that different behaviors in general involve different levels of arousal and movement, our model included behaviordependent modulation of frontal cortical neural activity.
Finally, for the reviewer’s convenience, we also quote below the paragraph addressing this issue in the revised Discussion.
“Another limitation of our model is the “openloop” nature of the relationship between behavior and neural activity. Specifically, we modeled neural activity as being modulated by behavior, but behavior was modeled using a Markov chain that is independent from the neural activity. In reality, neural activity and behavior form a closedloop, with different social behaviors being controlled by the neural activity of specific neural populations in specific brain regions. Thus, an important future direction is to close the loop by incorporating neural control of social behaviors into models of the interbrain relationship in bats. This will require future experimental studies to identify which frontal cortical regions and populations in bats are necessary or sufficient to control social behaviors, as well as the detailed causal relationship from neural activity to social behavior. Furthermore, as social interactions can occur at multiple timescales, it will be interesting to investigate how these are controlled by neural activity at different timescales, and how those timescales are shaped by functional acrossbrain coupling. In summary, such a closedloop model will shed light on how interbrain activity patterns and dynamic social interactions coevolve and feedback onto each other.”
Further, the Markov assumption is not rigorously tested.
We have now tested the Markov assumption, using the following methods. We compared three models of bat behaviors: (1) the independent model, where the behavioral state at a given time point is independent from the state at other time points; (2) the 1storder dependency model, where the behavioral state at a given time point depends on the state at the previous time point only; (3) the 2ndorder dependency model, where the behavioral state at a given time point depends on the states at the two previous time points. The Markov assumption corresponds to model (2), which is used as a part of the main model of the paper. Note that models with longer timedependencies (≥3) were not tested because the number of parameters grows exponentially with model order and our dataset is not large enough to fit them.
To compare the three models, we split the behavioral data into a training set and a test set, fitted each model on the training set (Laplace smoothing was used to avoid assigning zero probability to unobserved events), and calculated the loglikelihood of the test set under each model. Figure 3—figure supplement 2AB shows the crossvalidated likelihoods for the behavioral data of onechamber (A) and twochambers (B) sessions, which were fitted separately; circles and error bars are means and standard deviations across 100 random splits of the data into training and test sets.
As Figure 3—figure supplement 2AB shows, the 1storder model had the highest likelihood on average. This does not necessarily prove that bat behavior obeys the Markov assumption (if we had a lot more data, we might be able to fit better 2ndorder and higherorder models). But this does mean that, given the amount of data we have, the best model that we can fit is the 1storder Markov chain. Thus, this result supports our usage of the Markov chain in the main model of the paper.
In the revised manuscript, the analysis is described in Materials and methods.
No model selecting or other model validation appears to be done.
To evaluate model fit, we simulated our model using experimentally observed behaviors (rather than simulating behaviors using a Markov chain), and compared the simulated neural activity with the experimentally observed activity (see Materials and methods for detailed procedures). The comparison for an example experimental session is shown in Figure 3—figure supplement 1, where we’ve plotted the experimentally observed neural activity and behaviors for bat 1 (A) and bat 2 (B), along with the simulated neural activity. The correlation coefficient between data and model are indicated above each plot. These are representative examples, as the average correlation over all sessions and bats is 0.72 (standard deviation is 0.10).
In evaluating model fit, we realized that the model in the original manuscript produced outputs with a DC offset different from that of the data. Thus, in the revised manuscript, we added one more behavioral parameter ($b}_{constant$) that adjusts the DC offset, which is a parameter that reflects the effect of a baseline arousal level on neural activity (Materials and methods). Note that, since the only effect of this parameter is to adjust the DC offset of neural activity, it does not change any of the results in the paper.
In short, the model, while very interesting, is so complex that it is literally impossible to evaluate. The authors report literally no shortcomings of their model. They do not report parameter estimation methods. They do not report fitting errors or other model validation metrics. The only evaluation is whether it can produce certain outputs that are similar to biological data. While the latter is certainly important, all models are wrong, and it essential to have a model simple enough to understand, both in terms of how it works and how it fails.
The comments on the complexity of the model and on fitting errors have been addressed above. Regarding parameter estimation methods, they were described in Materials and methods, and we regret having neglected to directly reference it in the original manuscript. We now reference it in the legend of Figure 3A which is the first place to introduce the parameters. Briefly, the behavioral parameters ($b}_{resting$, $b}_{fighting$, etc.) were simply chosen to be the average neural activity during the respective behaviors from the data; the other parameters were chosen by hand to roughly match the levels of activity from the data, keeping within the parameter regime of $0<{\mathit{C}}_{I}<{\mathit{C}}_{S}$ identified from the analyses. As we showed above, these parameters provide a reasonable fit to the data.
The reason we chose the parameters heuristically in this way, rather than by minimizing some error objective, is the following. Our goal was to build a model that could qualitatively reproduce the experimental findings in a robust manner, that is, without finetuning of parameters. Thus, we analyzed the model to understand how model behaviors depend on the parameters, and to identify the parameter regime that reproduces the qualitative trends seen in the data (Figure 3IJ; Materials and methods). Guided by these analyses, we chose parameters heuristically without algorithmic finetuning.
Finally, following suggestions from reviewer 1 and reviewer 3, we have added discussions of shortcomings of the models (the last two paragraphs of the Discussion). With these discussions of model limitations, along with the presentation of simple insights into model mechanism from the reduced models above, we believe we have now presented a model that is “simple enough to understand, both in terms of how it works and how it fails.”
In general, while the basic finding is fairly interesting, and the experiments and their findings are highly relevant to the field, the modeling and its explication fall short.
It is not that it is wrong or bad; however, it is not clear that such a complex model increases our understanding beyond the experimental findings in Figure 1, and if it does, there has to be a major caveat that the model itself is not carefully vetted.
Based on the reviewer’s comments on the model’s complexity, we have analyzed reduced versions of the model to understand its simple underlying mechanisms, as described above. This goes beyond the experimental findings in Figure 1, as it provides a computational mechanism that could give rise to those experimental findings. Moreover, based on the reviewer’s comments, we have more carefully vetted the model, by evaluating model fit and testing different behavioral models that assume or doesn’t assume the Markov property. Finally, we now discuss caveats of the model in the Discussion section, including the openloop nature of the model as pointed out by the reviewer.
https://doi.org/10.7554/eLife.70493.sa2Article and author information
Author details
Funding
National Institutes of Health (DP2DC016163)
 Michael M Yartsev
National Institute of Mental Health (1R01MH2538701)
 Michael M Yartsev
New York Stem Cell Foundation (NYSCFRNI40)
 Michael M Yartsev
Alfred P. Sloan Foundation (FG20179646)
 Michael M Yartsev
Brain Research Foundation (BRFSG201709)
 Michael M Yartsev
National Science Foundation (NSF 1550818)
 Michael M Yartsev
David and Lucile Packard Foundation (201766825)
 Michael M Yartsev
Simons Foundation
 Michael M Yartsev
Pew Charitable Trusts (00029645)
 Michael M Yartsev
Dana Foundation
 Michael M Yartsev
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank ES Sevilla, A Raha, J Chau, K Moi, N Juthani, M Zuercher, L Kasraie, and C Tran for behavioral annotation; S A Afjei for histology; A Halley and L Krubitzer for providing the image of the 3Dreconstructed brain; L Jiang for machining; B Olshausen, W Liberti, and members of the Yartsev Lab for discussion and comments on the manuscript; C Ferrecchia and G Lawson for veterinary oversight; the staff of the Office of Laboratory Animal Care for support with animal husbandry and care. This research was supported by NIH (DP2DC016163), NIMH (1R01MH2538701), the New York Stem Cell Foundation (NYSCFRNI40), the Alfred P Sloan Foundation (FG2017–9646), the Brain Research Foundation (BRFSG2017–09), National Science Foundation (NSF 1550818), the Packard Fellowship (2017–66825), the KlingensteinSimons Fellowship, the Pew Charitable Trust (00029645), and the Dana Foundation (to MMY).
Ethics
All experimental procedures complied with all relevant ethical regulations for animal testing and research and were approved by the Institutional Animal Care and Use Committee of the University of California, Berkeley (protocol number AUP20150171222).
Senior Editor
 Ronald L Calabrese, Emory University, United States
Reviewing Editor
 Noah J Cowan, Johns Hopkins University, United States
Publication history
 Received: May 18, 2021
 Preprint posted: June 2, 2021 (view preprint)
 Accepted: February 8, 2022
 Accepted Manuscript published: February 10, 2022 (version 1)
 Version of Record published: March 24, 2022 (version 2)
Copyright
© 2022, Zhang et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 1,343
 Page views

 271
 Downloads

 1
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Genetics and Genomics
 Neuroscience
For at least two centuries, scientists have been enthralled by the “zombie” behaviors induced by mindcontrolling parasites. Despite this interest, the mechanistic bases of these uncanny processes have remained mostly a mystery. Here, we leverage the Entomophthora muscaeDrosophila melanogaster “zombie fly” system to reveal the mechanistic underpinnings of summit disease, a manipulated behavior evoked by many fungal parasites. Using a highthroughput approach to measure summiting, we discovered that summiting behavior is characterized by a burst of locomotion and requires the host circadian and neurosecretory systems, specifically DN1p circadian neurons, pars intercerebralis to corpora allata projecting (PICA) neurons and corpora allata (CA), the latter being solely responsible for juvenile hormone (JH) synthesis and release. Using a machine learning classifier to identify summiting animals in real time, we observed that PICA neurons and CA appeared intact in summiting animals, despite invasion of adjacent regions of the “zombie fly” brain by E. muscae cells and extensive host tissue damage in the body cavity. The bloodbrain barrier of flies late in their infection was significantly permeabilized, suggesting that factors in the hemolymph may have greater access to the central nervous system during summiting. Metabolomic analysis of hemolymph from summiting flies revealed differential abundance of several compounds compared to nonsummiting flies. Transfusing the hemolymph of summiting flies into nonsummiting recipients induced a burst of locomotion, demonstrating that factor(s) in the hemolymph likely cause summiting behavior. Altogether, our work reveals a neuromechanistic model for summiting wherein fungal cells perturb the fly’s hemolymph, activating a neurohormonal pathway linking clock neurons to juvenile hormone production in the CA, ultimately inducing locomotor activity in their host.

 Neuroscience
Deep brain stimulation targeting the posterior hypothalamus (pHypDBS) is being investigated as a treatment for refractory aggressive behavior, but its mechanisms of action remain elusive. We conducted an integrated imaging analysis of a large multicentre dataset, incorporating volume of activated tissue modeling, probabilistic mapping, normative connectomics, and atlasderived transcriptomics. Ninetyone percent of the patients responded positively to treatment, with a more striking improvement recorded in the pediatric population. Probabilistic mapping revealed an optimized surgical target within the posteriorinferiorlateral region of the posterior hypothalamic area. Normative connectomic analyses identified fiber tracts and functionally connected with brain areas associated with sensorimotor function, emotional regulation, and monoamine production. Functional connectivity between the target, periaqueductal gray and key limbic areas – together with patient age – were highly predictive of treatment outcome. Transcriptomic analysis showed that genes involved in mechanisms of aggressive behavior, neuronal communication, plasticity and neuroinflammation might underlie this functional network.