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

Experiments were conducted inside 40.6 × 33.7 × 52.1 cm (length × width × height) cages, placed in dark chambers. High-speed 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 one-chamber sessions and two-chambers sessions. In the one-chamber sessions, bats behaved freely and interacted with each other inside a chamber. In the two-chambers sessions, bats behaved freely in separate, identical chambers. There were three types of two-chambers 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.

For the two-bat 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 non-social. 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.

• Self-grooming. 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.

Anesthesia and surgery were performed as described previously in Yartsev and Ulanovsky, 2013 to implant the bats with four-tetrode 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.

Electrophysiology was performed using a wireless neural data logging system (Neurolog-16 and MouseLog-16, Deuteron Technologies). Tetrodes were advanced ventrally every one to two days (mostly by 20–160 µm), to record activity at different sites.

To determine the neural recording sites, we performed histology as previously described (Yartsev and Ulanovsky, 2013), using Nissl staining.

To obtain LFP, we low-pass filtered the raw voltage traces (cut-off frequency: 200 Hz), and then downsampled it to 496.6 Hz (the two-bat experiments) or 520.8 Hz (the four-bat experiments).

We detected spikes from band-pass 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.

To 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 peak-normalized the power at each frequency for each spectrogram.

Firing rates for single units and multiunits were computed in 5 s bins with 2.5 s overlap between consecutive bins.

Previously, 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 peak-normalized 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).

The statistical tests used are stated in the figure legends. A significance level of 0.05 was used for all tests. Tests were two-tailed unless otherwise indicated.

The neural activity of two bats can be represented as a two-dimensional vector ${\displaystyle \left(\begin{array}{c}{\mathit{a}}_1\left(t\right)\\ {\mathit{a}}_2\left(t\right)\end{array}\right)}$, where ${\displaystyle {\mathit{a}}_1\left(t\right)}$ and ${\displaystyle {\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: ${\displaystyle \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 ${\displaystyle {\mathit{a}}_M\left(t\right)}$ is the mean component of the activity, and ${\displaystyle {\mathit{a}}_D\left(t\right)}$ is the difference component. Here, the difference component is defined as ${\displaystyle \frac{1}{2}\left[{\mathit{a}}_1\left(t\right)-{\mathit{a}}_2\left(t\right)\right]}$ rather than ${\displaystyle {\mathit{a}}_1\left(t\right)-{\mathit{a}}_2\left(t\right)}$ so as to have the same scale as the mean component. ${\displaystyle {\mathit{a}}_M\left(t\right)}$ represents the common activity pattern shared between the brains, and ${\displaystyle {\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.

In our data from one-chamber 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 inter-brain 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 inter-brain correlations, large mean components, and small difference components?

To answer this, we examine the relationship between inter-brain correlation and the mean and difference components. Let’s use ${\displaystyle {\mathit{A}}_1}$, a column vector, to denote the activity of bat 1 during a session, and similarly, ${\displaystyle {\mathit{A}}_2}$ for the activity of bat 2. ${\displaystyle {\mathit{A}}_1}$ and ${\displaystyle {\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 ${\displaystyle {\mathit{A}}_M=\left({\mathit{A}}_1+{\mathit{A}}_2\right)/2}$, and the activity of the difference component is ${\displaystyle {\mathit{A}}_D=\left({\mathit{A}}_1-{\mathit{A}}_2\right)/2}$. We use ${\displaystyle {\hat{\mathit{A}}}_1}$ to denote the mean-subtracted version of ${\displaystyle {\mathit{A}}_1}$ (here ‘mean-subtracted’ 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, ${\displaystyle {\hat{\mathit{A}}}_2}$, ${\displaystyle {\hat{\mathit{A}}}_M}$, and ${\displaystyle {\hat{\mathit{A}}}_D}$ for mean-subtracted versions of the respective vectors. Then, the Pearson correlation coefficient between the activity of two brains is:

Thus, the correlation depends on the quantities ${\displaystyle {\hat{\mathit{A}}}_M^T{\hat{\mathit{A}}}_M}$, ${\displaystyle {\hat{\mathit{A}}}_D^T{\hat{\mathit{A}}}_D}$, and ${\displaystyle {\left({\hat{\mathit{A}}}_M^T{\hat{\mathit{A}}}_D\right)}^2}$. Note that ${\displaystyle {\hat{\mathit{A}}}_M^T{\hat{\mathit{A}}}_M}$ is $N$ times the variance of the mean component, and ${\displaystyle {\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 inter-brain 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 ${\displaystyle {\mathit{A}}_D}$, ${\displaystyle {\hat{\mathit{A}}}_D}$, etc. as above, and denote their surrogate counterparts by ${\displaystyle {\mathit{S}}_D}$, ${\displaystyle {\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 1000-second moving average filter, and subtracted from it the average across its elements. Let’s call the resulting vector ${\displaystyle \hat{\mathit{R}}}$. Due to the smoothing, ${\displaystyle \hat{\mathit{R}}}$ varied on slow timescales. We then found the component of ${\displaystyle \hat{\mathit{R}}}$ that is orthogonal to ${\displaystyle {\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 ${\displaystyle {\hat{\mathit{S}}}_D^{\prime }}$ as a linear combination of ${\displaystyle {\hat{\mathit{A}}}_M}$ and ${\displaystyle {\hat{\mathit{R}}}_O:{\hat{\mathit{S}}}_D^{\prime }=}$ ${\displaystyle \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}}$ ${\displaystyle +{\hat{\mathit{R}}}_O/\sqrt{{\hat{\mathit{R}}}_O^T{\hat{\mathit{R}}}_O}}$. The surrogate for ${\displaystyle {\hat{\mathit{A}}}_D}$ is the vector ${\displaystyle {\hat{\mathit{S}}}_D}$, obtained by scaling ${\displaystyle {\hat{\mathit{S}}}_D^{\prime }}$ to have the same vector norm as ${\displaystyle {\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 ${\displaystyle {\mathit{S}}_D={\hat{\mathit{S}}}_D+{\overline{\mathit{A}}}_D}$, where every element of the vector ${\displaystyle {\overline{\mathit{A}}}_D}$ is equal to the average across the elements of ${\displaystyle {\mathit{A}}_D}$. The surrogate activity of bat 1 and bat 2 are, respectively, ${\displaystyle {\mathit{S}}_1={\mathit{A}}_M+{\mathit{S}}_D}$ and ${\displaystyle {\mathit{S}}_2={\mathit{A}}_M-{\mathit{S}}_D}$. This procedure ensures that ${\displaystyle {\hat{\mathit{S}}}_D^T{\hat{\mathit{S}}}_D={\hat{\mathit{A}}}_D^T{\hat{\mathit{A}}}_D}$ and ${\displaystyle {\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 inter-brain 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 ${\displaystyle {\mathit{A}}_M}$ from the actual data unchanged in generating the surrogate data, so that the surrogate ${\displaystyle {\mathit{S}}_1}$ and ${\displaystyle {\mathit{S}}_2}$ show qualitatively similar behavioral modulation over time as the actual data ${\displaystyle {\mathit{A}}_1}$ and ${\displaystyle {\mathit{A}}_2}$, but it is also possible to replace both ${\displaystyle {\mathit{A}}_M}$ and ${\displaystyle {\mathit{A}}_D}$ with surrogate counterparts.

For Figure 2D–E, we repeated the above procedure for each of the one-chamber sessions from the actual data set, using the actual 30–150 Hz LFP power data.

Our 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; ${\displaystyle \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$; ${\displaystyle \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 ${\displaystyle {\mathit{C}}_S}$ is the strength of functional self-coupling and ${\displaystyle {\mathit{C}}_I}$ is the strength of functional across-brain coupling; ${\displaystyle \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 ${\displaystyle {\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 ${\displaystyle {\mathit{b}}_2\left(t\right)}$ .

The functional across-brain coupling term, ${\displaystyle {\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 one-chamber sessions, ${\displaystyle {\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 two-chambers sessions, ${\displaystyle {\mathit{C}}_I=0}$ since the two bats are separated. To ensure stability (so that neural activity do not go to infinity), ${\displaystyle -{\mathit{C}}_S}$ must be negative and must have a larger absolute value than ${\displaystyle {\mathit{C}}_I}$ ; thus, ${\displaystyle 0\le {\mathit{C}}_I<{\mathit{C}}_S}$ .

To generate ${\displaystyle \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 self-grooming would be one state. For the Markov chain for one-chamber 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 one-chamber 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 one-chamber 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 two-chambers Markov chain. The transition probability matrices for one-chamber and two-chambers sessions are shown in Figure 3—figure supplement 2C-D.

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, ${\displaystyle \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$, ${\displaystyle \mathit{b}\left(t_1\right)}$ was the sum of the noiseless deterministic behavioral modulation ${\displaystyle {\mathit{b}}_d\left(t_1\right)}$ and the behavioral modulation noise ${\displaystyle {\mathit{b}}_n\left(t_1\right)}$. ${\displaystyle {\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 self-grooming at time $t_1$, then ${\displaystyle {\mathit{b}}_d\left(t_1\right)=\left(\begin{array}{c}{b}_{resting}\\ b_{self-grooming}\end{array}\right)+b_{constant}}$, where $b_{resting}$ and $b_{self-grooming}$ are parameters specifying the level of behavioral modulation associated with the two behaviors, and $b_{constant}$ is a constant offset that differs between the one-chamber and two-chambers 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 ${\displaystyle {\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 ${\displaystyle \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, ${\displaystyle \mathit{b}\left(t\right)}$ at time points in-between were linearly interpolated.

Having simulated the behaviors and generated ${\displaystyle \mathit{b}\left(t\right)}$ for a 100 minute session, we set the initial condition ${\displaystyle \mathit{a}\left(0\right)}$ to be the fixed point under the initial noiseless behavioral modulation ${\displaystyle {\mathit{b}}_d\left(0\right)}$: ${\displaystyle \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 ${\displaystyle \mathit{a}\left(t\right)}$. The only differences between simulations of one-chamber and two-chambers sessions were the value of ${\displaystyle {\mathit{C}}_I}$ and the behavior transition probability matrix. For analyses and figures involving ${\displaystyle \mathit{a}\left(t\right)}$, ${\displaystyle \mathit{b}\left(t\right)}$, ${\displaystyle {\mathit{b}}_d\left(t\right)}$, and ${\displaystyle {\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).

As 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 first-order 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 second-order dependency model, where the behavioral state at a given time point depends on the states at the two previous time points. Models with longer time-dependencies (≥ 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 one-chamber and two-chambers 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 log-likelihood 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 log-likelihoods are shown in Figure 3—figure supplement 2A, B.

As shown in Figure 3—figure supplement 2A, B, the first-order 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 second-order 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 first-order Markov chain, supporting its use in our main model described in the previous section.

To 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 ${\displaystyle \mathit{b}\left(t\right)}$ from behaviors simulated using a Markov chain, we generated ${\displaystyle \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.

How 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 one-chamber model, the activity of bat 1 ${\displaystyle {\mathit{a}}_1\left(t\right)}$ influences the activity of bat 2 ${\displaystyle {\mathit{a}}_2\left(t\right)}$, which in turn feeds back on ${\displaystyle {\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 ${\displaystyle \mathit{C}}$. The eigenvectors of ${\displaystyle \mathit{C}}$ are ${\displaystyle \left(\begin{array}{c}1\\ 1\end{array}\right)}$ and ${\displaystyle \left(\begin{array}{c}1\\ -1\end{array}\right)}$, with respective eigenvalues ${\displaystyle -{\mathit{C}}_S+{\mathit{C}}_I}$ and ${\displaystyle -{\mathit{C}}_S-{\mathit{C}}_I}$. We define ${\displaystyle \mathit{V}}$ to be the matrix whose columns are the eigenvectors: ${\displaystyle \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: ${\displaystyle {\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: ${\displaystyle {\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 one-chamber model, ${\displaystyle {\mathit{C}}_I>0}$, so functional across-brain 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$, ${\displaystyle {\mathit{a}}_M\left(t\right)}$ is exponentially approaching ${\displaystyle \frac{{\mathit{b}}_M\left(t\right)}{{\mathit{C}}_S-{\mathit{C}}_I}}$ with effective time constant ${\displaystyle \frac{\tau }{{\mathit{C}}_S-{\mathit{C}}_I}}$, and ${\displaystyle {\mathit{a}}_D\left(t\right)}$ is exponentially approaching ${\displaystyle \frac{{\mathit{b}}_D\left(t\right)}{{\mathit{C}}_S+{\mathit{C}}_I}}$ with effective time constant ${\displaystyle \frac{\tau }{{\mathit{C}}_S+{\mathit{C}}_I}}$. In the one-chamber model, ${\displaystyle 0<{\mathit{C}}_I<{\mathit{C}}_S}$, so ${\displaystyle 0<{\mathit{C}}_S-{\mathit{C}}_I<{\mathit{C}}_S+{\mathit{C}}_I}$. This means that in the one-chamber model, relative to the two-chambers model where ${\displaystyle {\mathit{C}}_I=0}$, the positive feedback provided by functional across-brain 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 one-chamber and two-chambers models were the value of ${\displaystyle {\mathit{C}}_I}$ and the behavioral transition probability matrix. Having understood the effects of ${\displaystyle {\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 one-chamber and two-chambers 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 two-chambers 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 one-chamber model, same-behavior states (probability 0.58) are more likely than different-behavior states (Figure 3—figure supplement 2G). The noiseless behavioral modulation ${\displaystyle {\mathit{b}}_d\left(t\right)}$ depends on the behaviors of the two bats, while the behavioral modulation noise ${\displaystyle {\mathit{b}}_n\left(t\right)}$ does not. Through ${\displaystyle {\mathit{b}}_d\left(t\right)}$, same-behavior states contribute to behavioral modulation of the mean component, and the only noiseless behavioral modulation of the difference component comes from different-behavior states. In the one-chamber model, with increased coordinated behavioral modulation compared to the two-chambers 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 across-brain 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 across-brain 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 inter-brain 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 across-brain coupling is needed. To explicitly test this, we performed an additional set of simulations of the one-chamber model. Here, after we generated each instantiation of ${\displaystyle \mathit{b}\left(t\right)}$, we used it to simulate neural activity twice: once with functional across-brain coupling and once without. Note that the same ${\displaystyle \mathit{b}\left(t\right)}$ is used for both simulations, including the same ${\displaystyle {\mathit{b}}_d\left(t\right)}$ generated from the one-chamber Markov chain and the same instantiation of randomly generated ${\displaystyle {\mathit{b}}_n\left(t\right)}$. Figure 4C, D and G shows that the model with functional across-brain coupling reproduced the experimental observation: inter-brain correlation remained after removing periods of coordinated behaviors. On the other hand, without functional across-brain coupling, inter-brain correlation disappeared after removing coordinated behaviors (Figure 4E–G). Thus, we conclude that, in the framework of our model, opposite feedback by functional across-brain coupling is necessary to reproduce the data.

To 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 ${\displaystyle \mathit{b}\left(t\right)}$ based on simulated behavior, ${\displaystyle {\mathit{b}}_M\left(t\right)}$ and ${\displaystyle {\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 ${\displaystyle {\mathit{C}}_S}$ and ${\displaystyle {\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 ${\displaystyle {\lambda }_M=-{\mathit{C}}_S+{\mathit{C}}_I}$ and ${\displaystyle {\lambda }_D=-{\mathit{C}}_S-{\mathit{C}}_I}$ are the eigenvalues of ${\displaystyle \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, ${\tilde{f}}_{M,k}$ is the Fourier coefficient of $f_M\left(t\right)$ at frequency $k$: ${\displaystyle {\stackrel{\sim }{f}}_{M,k}=\frac{1}{\sqrt{T}}{\int }_0^{T}dt{f}_M\left(t\right)e^{-i\frac{2\pi kt}{T}}}$. The approximation of (10) and (11) assumes periodic convolutions and ${\displaystyle \mathit{a}\left(0\right)=0}$ in (6) and (7), but it is still accurate when ${\displaystyle \mathit{a}\left(0\right)}$ is smaller or comparable to the scale of ${\displaystyle \mathit{b}\left(t\right)}$ and $T$ is large compared to the effective time constants ${\displaystyle \frac{1}{{\mathit{C}}_S-{\mathit{C}}_I}}$ and ${\displaystyle \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, ${\displaystyle {\mathit{b}}_M}$ has a flat power spectrum: ${\displaystyle {\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|{\tilde{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 ${\displaystyle {\lambda }_M<0}$ and ${\displaystyle {\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 ${\displaystyle {\stackrel{\sim }{f}}_{M,k}=\frac{1}{\sqrt{T}}{\int }_0^{T}dt{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 ${\displaystyle {\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 ${\displaystyle {\lambda }_M<0}$ and ${\displaystyle {\lambda }_D<0}$, for large $T$, (17) can be approximated as

The results ${\displaystyle r_{Var}\approx \frac{{\mathit{C}}_S+{\mathit{C}}_I}{{\mathit{C}}_S-{\mathit{C}}_I}}$ and ${\displaystyle 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 ${\displaystyle 0<{\mathit{C}}_I<{\mathit{C}}_S}$. In this regime, consistent with the earlier analysis of the full model, ${\displaystyle {\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.

We explored an alternative mechanism for inter-brain 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 ${\displaystyle {\mathit{a}}_1\left(t\right)=\frac{sin\left({\theta }_1\left(t\right)\right)+1}{2}}$ and ${\displaystyle {\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: ${\displaystyle K>0}$ for simulations of one-chamber sessions, and $K=0$ for simulations of two-chambers sessions.

For each simulation, ${\omega }_1$ and ${\omega }_2$ were drawn from lognormal distributions with means ${\displaystyle {\overline{\omega }}_1}$ and ${\displaystyle {\overline{\omega }}_2}$, respectively, and the same standard deviation ${\sigma }_{\omega }$. ${\displaystyle {\overline{\omega }}_1}$, ${\displaystyle {\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 ${\displaystyle {\mathit{a}}_1\left(t\right)}$ and ${\displaystyle {\mathit{a}}_2\left(t\right)}$, which were plotted in Figure 4—figure supplement 1A-D and analyzed in Figure 4—figure supplement 1E-F. This shows that the phase-coupling mechanism of the Kuramoto model is able to reproduce inter-brain 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.

The neural dynamics in the model are governed by (2). To express (2) under a different basis, ${\displaystyle \mathit{a}\left(t\right)}$, ${\displaystyle \mathit{C}}$, and ${\displaystyle \mathit{b}\left(t\right)}$ are transformed to ${\displaystyle {\mathit{a}}^{\prime }\left(t\right)={\mathit{U}}^{-1}\mathit{a}\left(t\right)}$, ${\displaystyle {\mathit{C}}^{\prime }={\mathit{U}}^{-1}\mathit{C}\mathit{U}}$, and ${\displaystyle {\mathit{b}}^{\prime }\left(t\right)={\mathit{U}}^{-1}\mathit{b}\left(t\right)}$, respectively, where the columns of the matrix ${\displaystyle \mathit{U}}$ are the new basis vectors. If ${\displaystyle \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 counter-clockwise 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 ${\displaystyle {\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, self-grooming, 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.

Here 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 ${\displaystyle \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 ${\displaystyle \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 ${\displaystyle \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 ${\displaystyle {\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 ${\displaystyle \tau }$, ${\displaystyle {\left|{\stackrel{\sim }{\mathit{b}}}_{M,k}\right|}^2}$, and ${\displaystyle {\left|{\stackrel{\sim }{\mathit{b}}}_{D,k}\right|}^2}$. In this case, the neural decay functions are ${\displaystyle f_M\left(t\right)=e^{-\frac{{\mathit{C}}_S-{\mathit{C}}_I}{\tau }t}}$ and ${\displaystyle 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 ${\displaystyle {\mathit{B}}_{M,k}={\left|{\stackrel{\sim }{\mathit{b}}}_{M,k}\right|}^2}$ and $F_{M,k}={\left|{\tilde{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 ${\displaystyle k_t}$ such that ${\displaystyle \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 ${\displaystyle \mathit{B}}$ by adding ${\displaystyle \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 ${\displaystyle \delta \mathit{B}}$.

Considering $r_{Var}$ and $r_{PSC}$ as functions of ${\displaystyle \mathit{B}}$, their gradients are