Interactions across hemispheres in prefrontal cortex reflect global cognitive processing

Domain: neuroscience
Published: 2026-04-11
Editorial outcome: minor_revision

Abstract
Brain functions involve processing in local networks as well as modulation from brain-wide signals, such as arousal. Dissecting the contributions of populations of neurons to these functions requires knowledge of interactions between brain areas. We investigated these interactions using dual hemisphere recordings of prefrontal cortex in monkeys performing a spatial memory task. To tease apart global processing from local interactions, we applied a novel statistical approach called pCCA-FA (a combination of probabilistic canonical correlation analysis and factor analysis) to analyze trial-to-trial variability in neuronal activity. We found substantial shared variability among neurons within each population, much of which was actually shared across populations and linked to an arousal process. Our work presents a path by which we can leverage multi-area recordings to reveal aspects of brain functions that are hidden in single-area recordings.

Introduction
Brain functions such as attention [1, 2], decision making [3], and working memory [4, 5] rely on coordination among neurons within brain areas as well as across brain areas. Thus, the spiking activity of individual neurons can be influenced by neurons from the same brain area, neurons in other brain areas, global modulatory signals, and more [6]. Teasing apart the impact of each of these influences on neural activity would enable a deeper understanding of how neurons in different brain areas interact to give rise to such brain functions.

To dissect the various influences on neural activity, we can consider how the activity of one neuron covaries with that of other nearby neurons. Many studies have examined how the activity of pairs of neurons in a single brain area covaries from trial-to-trial by computing spike count correlation (rsc) [7, 8]. These pairwise correlations have been widely used as a window to investigate a range of phenomena from cognitive states to circuit architectures [9, 10]. To go beyond pairs of neurons, a related body of work has identified the population-wide covariability structure in single brain areas. These population-level approaches have been used to study the trial-to-trial variability underlying sensory encoding [11, 12], attention [13, 14], motor control [15–17], learning [18–20], decision-making [21–23], network structure [24–26], and more [27]. Taken together, these studies inform our understanding of how neural activity is shared among neurons within a brain area.

Given that most brain functions rely on interactions across brain areas, it is likely that some component of the activity of each brain area is shaped by interactions with other brain areas. Developments in neural recording technology have enabled simultaneous recordings of populations of neurons in multiple brain areas [28–30], which give rise to the possibility of distinguishing interactions across brain areas from interactions solely within a single brain area [31–34]. Simultaneous recordings from multiple brain areas have been leveraged to study widespread brain functions such as sensory processing [35–40], motor control [41–44], learning [45–47], decision-making [50–52], working memory [48, 49], arousal [53, 54], and attention [55–60]. However, most statistical methods for analyzing populations of neurons are not designed to distinguish activity shared across areas from activity shared solely among neurons within each area.

Here, we seek to identify components of neural activity that involve interactions across two areas (i.e., an across-area component), as well as interactions solely within each area (i.e., a within-area component). We simultaneously recorded spiking activity from each of the two hemispheres (here referred to as two areas) of prefrontal cortex (PFC) in macaque monkeys while they performed a spatial memory task. Utilizing a statistical method we developed called pCCA-FA (a combination of probabilistic canonical correlation analysis and factor analysis), we found substantial shared variance among neurons in PFC, much of which involved neurons across both areas. In fact, the across-area component of PFC activity was larger than the within-area component. We then uncovered a link between the across-area component and a brainwide arousal process. Our work motivates a reinterpretation of shared variance among neurons within a brain area as arising from multiple sources, both within that area as well as via interactions across the brain.

Results
To support perception, cognition, and action, each neuron’s activity can be shared with other neurons in the same brain area, as well as with neurons in other brain areas. Consider the spiking activity of two distinct populations of neurons, recorded simultaneously from two brain areas (area 1 and area 2, Fig. 1a). For each neuron, we would like to identify which aspects of its activity are shared with neurons in the other area (pink arrow) and which aspects are shared exclusively with neurons in the same area (blue arrow). To understand how the activity of neuron 1 in area 1 is shared with other neurons (i.e., to quantify the relationships described by the pink and blue arrows in Fig. 1a), we would like to partition its trial-to-trial spike count variance into an across-area component (termed across-area variance; Fig. 1b, pink), a within-area component (termed within-area variance; Fig. 1b, blue), and a component independent to this neuron (termed independent variance; Fig. 1b, gray).

2.1 Mean spike count correlation does not adequately partition shared variance

We simultaneously recorded spiking activity from neurons in both hemispheres of dorsolateral PFC to understand how neural activity was shared either across or within areas (Fig. 2a). We began by characterizing the pairwise interactions among neurons within each area using rsc. The spike counts used to compute rsc were taken from a one-second window during the delay period of a memory-guided saccade task (Fig. 2b). Throughout this work, we combined trials across conditions after removing condition-specific means, and removed slow-timescale trends to avoid spurious correlations [61] (see Methods, Supplementary Fig. 1).

The mean rsc among pairs of recorded neurons in a brain area is commonly taken to indicate the strength of shared activity among neurons. We computed rsc for each pair of neurons within each area (Fig. 2c, left and middle, example session), then averaged over the pairs to obtain a within-area mean rsc per session (right PFC within-area mean rsc = 0.029, left PFC within-area mean rsc = 0.033, example session). In the majority of sessions, we found that within-area mean rsc values were positive (Fig. 2d, top and middle).

To understand how activity is shared across areas, one approach is to compute rsc between pairs of neurons across areas (i.e., one neuron in each area). Mean rsc decreases for pairs of neurons as the physical distance between the neurons increases [62–64], and has been shown to be smaller for pairs of neurons across versus within brain areas [65]. We observed an across-area mean rsc that was close to zero (Fig. 2c, right; Fig. 2d, bottom). This might be taken to suggest that neurons within areas interact more strongly than across areas (Fig. 2d; p = 1.0e−13, pooled across animals and brain areas; one-tailed t-test), and that there are little to no interactions across the hemispheres of PFC.

However, there are several limitations to mean rsc measurements. First, by averaging over all pairs of neurons (either across or within areas), we ignore the breadth of rsc distributions (Fig. 2c). For example, if there are pairs of neurons with strong positive correlations and also pairs with strong negative correlations, the mean rsc might be close to zero even though many neuron pairs exhibit strong correlations. One way to identify neuron pairs with strong positive and negative correlations is to consider the relative tuning preferences of the neurons [66–70]. Indeed, when we computed mean rsc among subsets of across-area neuron pairs with similar tuning (Fig. 2e, top) and opposite tuning (Fig. 2e, bottom), we found that similarly tuned pairs of neurons exhibited positive rsc (mean rsc = 0.019; p = 1.6e−10, pooled across animals; one-tailed t-test), whereas oppositely tuned pairs of neurons exhibited negative rsc (mean rsc = −0.013; p = 1.e−7, pooled across animals; one-tailed t-test). Second, the within-area mean rsc does not reveal whether the activity shared among neurons in an area is also shared with other areas. Rather, all co-fluctuations among neurons in a brain area are collapsed into the within-area mean rsc. As such, within-area mean rsc could reflect interactions with other areas and/or global neuromodulation. Third, mean rsc (either across- or within-area) is not designed to identify the many possible population-wide co-fluctuations among neurons. Several distinct changes in population-wide co-fluctuations can result in the same change in mean rsc [71, 72]. By collapsing all interactions in one metric, we are unable to separately characterize the many influences on neural activity. Here we desire to leverage the simultaneity of our dual PFC recordings to dissect the multiple population-wide influences on across- and within-area activity (as in Fig. 1b).

2.2 pCCA-FA partitions across- and within-area shared variance

To explicitly partition across- and within-area variance, we can leverage statistical methods that take into account all simultaneously-recorded neurons together. In a single brain area, we can identify population-wide co-fluctuation patterns using factor analysis (FA) or related methods [73]. While this provides a more detailed view of population-wide covariability than mean rsc [72], these methods do not reveal whether activity fluctuations are shared with other brain areas. When recording from two brain areas simultaneously, we can identify variance shared among neurons across areas by employing canonical correlation analysis (CCA) or related methods [31]. However, neither CCA nor its probabilistic variant, pCCA [74], are designed to also identify variance shared only among neurons within each area (see Methods). Given that none of the aforementioned methods are designed to partition variance into both across-area and within-area components, a new approach is needed.

To address this need, we developed a dimensionality reduction method for two simultaneously recorded brain areas called pCCA-FA (a combination of probabilistic canonical correlation analysis and factor analysis), which partitions the variance of each area’s activity into an across-area component, a within-area component, and a component independent to each neuron. pCCA-FA can be applied to any two populations of simultaneously recorded neurons, and here we focus on the case where each population is recorded in a distinct brain area (Fig. 3a, area 1 and area 2). We seek to identify latent variables that capture co-fluctuations among neurons across areas (Fig. 3a, z, pink), as well as latent variables specific to each area that capture co-fluctuations exclusively among neurons within each area (Fig. 3a, area 1: z1, light blue; area 2: z2, dark blue). In total, pCCA-FA identifies three sets of latent variables: one across-area set, and a separate within-area set for each area. Because the within-area latent variables are specific to each area, the number of within-area latent variables for each area need not be the same.

pCCA-FA approximates neural activity as a weighted sum of the across-area latent variables (Fig. 3b, pink) and the within-area latent variables (Fig. 3b, light and dark blue for areas 1 and 2, respectively). The co-fluctuation patterns comprise the weights, which are estimated from the neural activity and fixed across all trials. The values of the latent variables can vary from trial to trial. Each co-fluctuation pattern can be represented as a vector in a population activity space, where each coordinate axis represents the activity of one neuron (Fig. 3c). Whereas the across-area co-fluctuation patterns are paired across areas (Fig. 3c, pink), the within-area co-fluctuation patterns are separate for each area (i.e., not paired across areas) (Fig. 3c, light and dark blue for areas 1 and 2, respectively).

For each neuron, we can quantify the proportion of its spike count variance explained by each set of latent variables. pCCA-FA partitions each neuron’s variance into a component that is described by the across-area latent variables (Fig. 3d, pink) and a component that is described by that area’s within-area latent variables (Fig. 3d, light and dark blue for areas 1 and 2, respectively). The residual variance not accounted for by across- or within-area latent variables is independent to each neuron (Fig. 3d, gray). Thus, pCCA-FA explicitly partitions the spike count variance of each neuron into across-area, within-area, and independent components (as in Fig. 1b).

To summarize the outputs of pCCA-FA, we define the following two metrics. First, to understand the extent to which a particular set of latent variables explains a neuron’s spike count variance, we define the across- and within-area percentage of variance shared among neurons (%sv). The values of this metric reflect the relative sizes of the bars in Fig. 3d. Across-area %sv is defined as the proportion of spike count variance explained by the across-area latent variables, averaged across neurons in each brain area (see Equation (8)). Within-area %sv is defined analogously using the within-area latent variables (see Equation (9)). Thus, each brain area has both an across-area %sv and a within-area %sv. Second, we quantify the complexity of the interactions described by each set of latent variables by computing its shared dimensionality (dshared) [25, 75]. A larger dshared would indicate more complex interactions among neurons, whereas a smaller dshared would indicate less complex interactions among neurons. Across- and within-area dshared are defined as the number of latent variables required to explain the across- or within-area variance, respectively (see Methods).

Before applying pCCA-FA to neural recordings, we performed extensive model validation. We verified that pCCA-FA accurately identified the ground truth %sv and dshared using simulated data with properties chosen to be similar to our neural recordings (Supplementary Figs. 2-3). We also found that pCCA-FA required fewer trials to recover the ground truth %sv and dshared than pCCA (Fig. 3e-f) because pCCA-FA has fewer parameters than pCCA (Supplementary Fig. 4). Further, pCCA-FA outperformed other dimensionality reduction techniques (pCCA and FA) on neural recordings (Supplementary Fig. 4). Thus, the pCCA-FA framework successfully captures how neurons in two brain areas covary by partitioning variance into across- and within-area components.

2.3 Across-area interactions were substantial and often greater than within-area interactions

To identify across- and within-area variance in our PFC recordings, we fit pCCA-FA to the same neural activity (with the same preprocessing) as in the rsc analyses (Fig. 2). For each session, pCCA-FA estimated a single set of across-area latent variables common to both areas and a separate set of within-area latent variables for each area. From the mean rsc values in Fig. 2, one might expect to uncover a large within-area component and small across-area component using pCCA-FA.

To determine if this was the case, we computed across- and within-area %sv and dshared. For the same example session as in Fig. 2c, we found that a substantial portion of variance was shared across areas (14.97% and 17.47%, for right and left PFC, respectively), compared to the variance shared exclusively within each area (5.65% and 7.47%, for right and left PFC, respectively) (Fig. 4a). This trend was evident across sessions, with the across-area %sv being larger than the within-area %sv (Fig. 4b; p = 2.6e−7, pooled across animals and brain areas; paired one-tailed t-test; see Supplementary Fig. 5 for a shuffle control). Thus, substantial across-area interactions were present in the PFC activity and identified by pCCA-FA (see also Supplementary Fig. 1e for the across-area correlation values computed using pCCA-FA), contrary to the intuition from measuring pairwise across-area rsc. This discrepancy can be understood by considering the sign of pairwise correlations among neurons (Fig. 2e, see Math Note A of ref. [72]). Mean rsc averages over positive and negative correlations found in individual pairs of neurons, whereas %sv accounts for both positive and negative correlations among neurons. Note that the values of across- and within-area %sv (like mean rsc) can depend on neuronal properties, including their physical distance and tuning preferences.

Having found that neurons tended to more strongly covary across areas than within each area, we next considered the complexity of these interactions. We found that the across-area dshared was also larger than the within-area dshared (Fig. 4c; p = 3.1e−6, pooled across animals and brain areas; paired one-tailed t-test). This indicated that, in addition to being stronger, the interactions we identified across areas were more complex than the interactions within each area. We further note that these findings are not a consequence of statistical bias in measuring across- versus within-area %sv and dshared (Supplementary Fig. 2).

2.4 Across- and within-area components involved different co-fluctuation patterns

We next sought to understand the population-wide co-fluctuation patterns represented by each set of latent variables, as they provide insight into the nature of the co-fluctuations among neurons (Fig. 5a). There are many possible population-wide co-fluctuations. Perhaps the simplest example is all neurons sharing the same-signed loading onto a latent variable, such that the population would all increase their activity together or decrease their activity together. We examined whether this was the case in our PFC recordings. For this analysis, we focused on the first co-fluctuation pattern of each ordered set, as it reflected a large proportion of the across- or within-area variance (Supplementary Fig. 6). We found that all neurons increasing their activity together or decreasing their activity together was not consistent with the across- or within-area co-fluctuation patterns we identified (Fig. 5b, Supplementary Fig. 7).

We next asked whether across- and within-area co-fluctuation patterns were similar to each other by computing the angle between them in population activity space (θestimated; Fig. 5a). They might be similar if two different processes (one global and one local) had a similar effect among the neurons. Alternatively, each process might engage distinct circuit mechanisms, yielding different co-fluctuation patterns. In an example session, we found that the across- and within-area co-fluctuation patterns appeared visually distinct, and produced a θestimated of 72◦ (Fig. 5b). In all sessions, we found that the across- and within-area co-fluctuation patterns had large θestimated (Fig. 5c), often close to 90◦, which would indicate distinct co-fluctuation patterns.

To further assess the extent to which the measured θestimated reflected true differences in the across- and within-area co-fluctuation patterns, we created a simulation in which we defined the ground truth angle between across- and within-area co-fluctuation patterns (θsim, see Methods, Supplementary Fig. 3). We simulated activity in each area with a specified value of θsim that was otherwise matched to the neural recordings in terms of the number of neurons, latent variables, and trials (see Methods). We found that, for any given θsim, the resulting distributions of θestimated (Fig. 5d, top) were narrower than what we observed in the neural recordings (Fig. 5c). Instead, if we sampled a different value of θsim for each simulated dataset uniformly from 45◦ to 90◦, we recovered a distribution of θestimated (Fig. 5d, bottom) that was similar to that observed in the neural recordings (p = 0.40, pooled across animals and brain areas; paired two-tailed t-test). We confirmed that the observed variation in θestimated from session to session was not related to the number of neurons or trials recorded (Supplementary Fig. 8). Furthermore, we determined that the variation in θestimated did not arise solely from sampling different subsets of neurons from session to session, although it still likely contributed (Supplementary Fig. 9). Together, these results indicate that across- and within-area components involved different co-fluctuation patterns that varied in their similarity from session to session.

2.5 Across-area latent variables were related to pupil diameter

Given that we found substantial across-area variance (Fig. 4) involving different co-fluctuation patterns than the within-area component (Fig. 5), we considered whether the across-area variance might be related to a brainwide process. Arousal is known to modulate activity of neurons throughout the brain [53, 54, 76–78]. In addition, arousal state has been linked to changes in both the tonic diameter of the pupil [79, 80] and event-related changes in pupil diameter [81, 82]. Changes in pupil diameter have also been linked to arousal-related changes in neural activity in both cortical [83–86] and subcortical [79, 87–89] regions. Given that arousal has been found to be a brainwide phenomenon, we wondered if it was related to the across-area variance we observed. Recall that the across-area component is defined as variance that is shared with one or more neurons in the other area (regardless of whether that variance is also shared with neurons in the same area). We therefore asked how well the across- and within-area latent variables identified by pCCA-FA explained variance in the pupil diameter (Fig. 6a).

To address this question, we measured the trial-by-trial pupil diameter during the delay period and adjusted the pupil diameter values to remove slow fluctuations in an analogous manner to the neural activity, yielding one pupil diameter measure per trial (see Methods, Supplementary Fig. 10a-b). We used the across- or within-area latent variables identified by pCCA-FA (i.e., either d across-area latent variables or dm within-area latent variables for area m; see Methods) to compute a trial-by-trial regression of the pupil diameter for each set of latent variables (Fig. 6b). Across sessions, we found that the across-area latent variables explained substantially more variance in the pupil diameter (mean r2 = 0.25) than the within-area latent variables (right PFC, mean r2 = 0.038; left PFC mean r2 = 0.030) (Fig. 6c; p = 1.3e−30, pooled across animals and brain areas; one-tailed t-test). Furthermore, we found that the across-area latent variables explained variance in the event-related changes in pupil diameter observed just after target onset (Supplementary Fig. 10; p = 1.1e−5, pooled across animals and brain areas; one-tailed t-test). Thus, the across-area latent variables may be related to multiple aspects of pupil-linked arousal. Taken together, these results are consistent with a link between the coordinated activity across brain areas and global cognitive processes such as arousal.

Methods
4.1 Animal model
Three male rhesus macaque monkeys (Macaca mulatta) were used in this study. One animal was sourced from Tulane National Primate Research Center (Monkey Wa) and two animals were from Covance (Monkey Pe and Monkey Sa). They were 10 (Monkey Pe), 7 (Monkey Wa), and 7 (Monkey Sa) years old at the time of data collection. All experimental procedures were approved by the Institutional Animal Care and Use Committees of Carnegie Mellon University and the University of Pittsburgh, and were in accordance with the United States National Research Council’s Guide for the Care and Use of Laboratory Animals.

4.2 Surgical preparation
Each animal was surgically implanted with a titanium headpost, which was fixed onto the skull using titanium screws. This was necessary to limit head movement such that the position of the eyes could be tracked and neural activity recorded. Animals were then each implanted with two 96-electrode “Utah” arrays (Blackrock Neurotech, Salt Lake City, UT). Electrode arrays were placed on the prearcuate gyrus just anterior to the arcuate sulcus, one in each hemisphere. Surgeries were performed in sterile conditions under general anesthesia using isoflurane.

4.3 Electrophysiological methods
Animals were positioned 36 cm from a 21” cathode ray tube monitor with 1024x768-pixel resolution and a refresh rate of 100 Hz. Task stimuli were generated using custom-designed software in Matlab (Mathworks, Natick, MA) that utilized the Psychophysics toolbox extensions [139–141]. Signals from the implanted electrodes were band-pass filtered (0.3 - 7500 Hz) and then digitized at 30,000 Hz. For each electrode, spiking waveforms were defined as a 52-sample (1.73 ms) window of the filtered voltage signal triggered by the signal crossing a predefined threshold, and stored for offline analysis. The threshold was defined as a multiple of the root-mean-square voltage of a brief epoch of the raw signal on each electrode collected at the beginning of each session. All behavioral data and neural activity were recorded using a Grapevine recording system (Ripple, Salt Lake City, UT) for further offline processing and analysis.

4.4 Behavioral task
Animals were trained to perform a standard memory-guided saccade task [142]. At the beginning of each trial, a 0.5 degree diameter blue fixation circle appeared at the center of a gray screen. The animal initiated fixation within an invisible window (Monkey Pe: 2.4, Monkey Wa: 1.8, Monkey Sa: 2.9 degrees radius) centered on the fixation circle and then 200 ms later a white target flashed in the animal’s periphery (Monkey Pe: 11.8 or 17.5 degrees from fixation depending on the session, Monkey Wa: 11.6 degrees from fixation, Monkey Sa: a single value in the range of 8.0 - 15.6 degrees from fixation depending on the session, with one session containing multiple randomly interleaved distances in this range) at one of 4, 8, or 16 locations depending on the session. The white target remained on the screen for either 100, 200, or 400 ms (depending on the session) after which it was removed from the screen. The animal maintained fixation for a delay period drawn randomly from a uniform distribution between 1.5 and 3 seconds. After the delay, the fixation circle disappeared, which served as the go cue for the animal to initiate a saccade toward the remembered location where the white target flash had appeared. The animal had 500 ms (Monkey Sa: 400 ms) to initiate a saccade, and an additional 200 ms to reach the target location. Successful completion of a saccade was defined by the animal maintaining fixation within an invisible window (Monkey Pe: 4.1, Monkey Wa: 5.9 or 7.0, Monkey Sa: 4.1 degrees radius) centered at the target location for at least 100 ms. A liquid reward was delivered for a saccade to the correct location. For a subset of sessions, a dim white circle was flashed at the target location after saccade initiation to assist the animal in target acquisition. Trials were pseudo-randomized in blocks during which the animal was required to correctly complete all target directions before beginning a new block. In total, we collected 42 sessions of data (16 in Monkey Pe, 10 in Monkey Wa, and 16 in Monkey Sa). For each session, we considered only trials in which the animal successfully acquired the saccade target (i.e., correct trials). The number of correct trials per session for Monkey Pe was 729.8 ± 173.4 (mean ± standard deviation), for Monkey Wa was 463.5 ± 170.1, and for Monkey Sa was 850.1 ± 284.2.

4.5 Preprocessing of neural activity
We used a previously trained and published neural network to classify waveforms as spikes or noise [143]. Briefly, the network was trained on human-sorted waveforms to distinguish between waveforms putatively of neural origin and waveforms not of neural origin. Using this network, we removed threshold crossings that were unlikely to be of neural origin, and refer to the remaining waveforms as spikes. The spikes from each electrode on the Utah array were treated as multi-unit activity, and no further differentiation of single units was performed. Thus, each electrode was treated as one neural unit. Using multi-unit activity tends to yield higher correlation values than single unit activity [8], and this applies equally to the across-area and within-area measurements.

We removed electrodes that had very low firing rates or very high trial-to-trial variability. To do this, we first binned neural activity by counting spikes during the target presentation and delay period between target onset and saccade initiation cue. We then removed electrodes with mean spike count lower than 2 spikes/second or Fano factor greater than 10 during this time window. This served to eliminate electrodes contaminated by recording artifacts, or units that were not present for the entirety of a recording session. We also removed electrodes affected by artifactual cross-talk due to electrical coupling. For each pair of electrodes, we flagged spikes as coincident if they occurred within 100 µs of each other. If either electrode in the pair had 20% or more of its spikes flagged as coincident, we flagged that pair as having artifactual crosstalk. This procedure provided an index for each electrode of how many other electrodes it was coincident with. We then iteratively removed the electrode coincident with the highest number of other electrodes until no crosstalk remained, which ensured that the fewest possible number of electrodes were removed. After this process the number of remaining electrodes in Monkey Pe was 78.3 ± 7.8 and 79.5 ± 8.5 (mean ± standard deviation), in Monkey Wa was 85.3 ± 8.1 and 25.2 ± 4.4, and in Monkey Sa was 62.6 ± 9.8 and 75.3 ± 13.2 for right and left PFC respectively.

After the above preprocessing steps, the spiking activity was binned into spike counts. The bin was defined as the one second of the delay period immediately preceding the earliest possible saccade initiation cue in the session, such that each trial contributed one bin. For example, if the earliest possible saccade initiation cue in a given session was 1.5s after the target disappeared, the bin would be defined as 0.5-1.5s after the target disappeared for every trial in that session. This one-second window aggregated the interactions across and within brain areas observed during the delay period on the trial-to-trial timescale.

4.6 Computing signal correlation
In Fig. 2e, we investigated if the distribution of across-area rsc, despite having a mean near zero, contained pairs of neurons with positive rsc and other pairs of neurons with negative rsc, in a way that was systematically related to the neurons’ task tuning. To do this, we identified each neuron’s tuning across all of the saccade target locations. We determined whether each pair of neurons had similar or opposite tuning by computing their signal correlation. This was defined as the Pearson correlation between the two neurons’ raw trial-averaged responses to each of the possible targets (4, 8, or 16 depending on the session) during the one-second window described above. To assess significance for each neuron pair, we generated a chance distribution of correlation values by using a permutation test (shuffling the target location labels). We then labeled a neuron pair as having significant signal correlation (positive for similarly tuned and negative for oppositely tuned pairs) if their signal correlation was more extreme than the 99th percentile of the chance distribution. Neuron pairs not meeting this criterion were excluded from this analysis.

4.7 Combining trials across conditions
To increase the number of observations for model fitting, we combined trials across conditions (saccade target direction). To focus our analyses on trial-to-trial variability, for each neuron, we first subtracted the condition mean from the binned spike counts in each condition and then combined residual spike counts across conditions. Because the spike count variance can differ across conditions, the residual spike counts can still contain a small amount of condition-specific signal. An alternative approach is to z-score the spike counts for each neuron and condition before combining across conditions [63]. We repeated the key analyses related to the identification of shared variance among neurons in Figs. 2 and 4 using condition-specific z-scoring and found qualitatively similar results. Note, the computation of signal correlation (above) was performed before the subtraction of condition-specific means.

4.8 Separation of slow and fast components
In this work, we were interested in trial-to-trial variability of neural activity on a timescale of seconds. As we investigated interactions between left and right PFC, it soon became apparent that both areas also contained a component that varied slowly on a timescale of minutes to hours. While slow-timescale variability presents its own set of scientific questions (e.g., ref. [54]), these fluctuations complicated our analyses by introducing dependencies between one trial and the next, violating the basic independence assumption of correlation analysis (including regression, Pearson correlation, pCCA, and pCCA-FA). Therefore, we removed the slow component from the neural activity (see Supplementary Fig. 1), and focused all analyses in this work on faster-timescale trial-to-trial variability, except where noted in Supplementary Fig. 1 and for computing signal correlation. We computed the slow component for each session using a centered boxcar filter of length 25 trials, computed for each remaining electrode after the preprocessing described above. We then subtracted this component from the residual spike counts to remove slow-timescale trends that could have induced spurious correlations [61]. We performed this preprocessing procedure independently for each neuron. All analyses in this study were performed on the residual faster-timescale component, except where noted.

4.9 Preprocessing of pupil diameter
Eye position and pupil diameter were monitored using monocular infrared tracking at a 1000 Hz sample rate (EyeLink 1000, SR Research, Mississauga, Canada). On each trial, we measured the pupil diameter during the delay period, from target offset until the saccade initiation cue. Each trial was smoothed using a median filter with a 50ms window (i.e., the median value in a 50ms centered window was subtracted from each sample) to remove high frequency noise.

For the analyses in Fig. 6, we computed mean smoothed pupil diameter (arbitrary units) over the full time window described above to obtain the most robust estimate of the pupil diameter on each trial. This yielded one average pupil diameter value per trial. Like the preprocessing of neural activity, we removed slow-timescale fluctuations in pupil diameter using a centered boxcar filter of length 25 trials. This yielded what we termed the trial-by-trial pupil diameter (Supplementary Fig. 10a-b).

To ensure pupil data quality, we computed the variance across trials of the trial-by-trial pupil diameter. We then removed outlier sessions that exceeded a variance threshold of two times the interquartile range of the variance values across sessions. This was performed for each animal individually. This procedure excluded one session from Monkey Wa from analyses of the pupil (Fig. 6, Supplementary Fig. 10), resulting in n = 9 sessions for analysis. To ensure pupil diameter was on a similar scale across sessions, we z-scored the trial-by-trial pupil diameter for each session by subtracting the mean and dividing by the standard deviation across trials in the session.

4.10 Probabilistic canonical correlation analysis (pCCA)
A commonly used approach for identifying interactions between brain areas is canonical correlation analysis (CCA) [31]. Briefly, CCA attempts to find a co-fluctuation pattern in area 1 and a co-fluctuation pattern in area 2, such that neural activity along these co-fluctuation patterns is maximally correlated across areas. Further co-fluctuation patterns can be found subject to the constraint that they are uncorrelated with previous co-fluctuation patterns. The correlations along these co-fluctuation patterns are known as canonical correlations (ρ). Note that CCA identifies activity shared across areas, but is not designed to partition activity into across- and within-area components.

The probabilistic interpretation of CCA (pCCA [74]) is defined by the following model. pCCA defines a linear-Gaussian relationship between the spike count vectors x1 ∈ Rn1 from area 1 and x2 ∈ Rn2 from area 2 and the across-area latent variables z ∈ Rd as:

z ∼ N(0, Id)
x1 | z ∼ N(µ1 + W1z, Ψ1)
x2 | z ∼ N(µ2 + W2z, Ψ2) (1)

where n1 and n2 are the number of neurons recorded in areas 1 and 2, respectively, and d is the number of across-area latent variables. The identity matrix Id ∈ Rd×d describes prior covariance of the latent variables. The vectors µ1 ∈ Rn1 and µ2 ∈ Rn2 represent the mean spike counts of areas 1 and 2, respectively. The loading matrices W1 ∈ Rn1×d and W2 ∈ Rn2×d map the latent variables onto neural activity in areas 1 and 2 respectively, and Ψ1 ∈ Rn1×n1 and Ψ2 ∈ Rn2×n2 are full covariance matrices.

From Equation (1), the marginal distributions for x1 and x2 are:

x1 ∼ N(µ1, W1W1T + Ψ1)
x2 ∼ N(µ2, W2W2T + Ψ2) (2)

We can see that pCCA decomposes the covariance of neural activity from each area into a low-rank component shared across areas (e.g., W1W1T) and a full noise covariance matrix (e.g., Ψ1). Note that pCCA treats all variance which cannot be explained by across-area latent variables (z), including within-area interactions and variance independent to each neuron, as noise.

It can be shown that pCCA (Equation (1)) finds the same co-fluctuation patterns (described by W1 and W2) as CCA that maximize correlation across areas [74]. CCA and pCCA therefore identify the same canonical correlations.

4.11 Probabilistic canonical correlation analysis - factor analysis (pCCA-FA)
We developed a model called pCCA-FA (a combination of probabilistic canonical correlation analysis and factor analysis (FA [144])), which explicitly partitions covariance of neural activity into an across-area component, a within-area component, and a component independent to each neuron. Though we term the components across-area and within-area, the pCCA-FA framework can be applied to any two simultaneously-recorded populations of neurons. pCCA-FA is related to extensions of FA that have been proposed in other contexts, such as refs. [145] and [146].

pCCA-FA defines the linear-Gaussian relationship between the spike count vectors x1 ∈ Rn1 from area 1 and x2 ∈ Rn2 from area 2 and across-area latent variables z ∈ Rd and within-area latent variables z1 ∈ Rd1 and z2 ∈ Rd2 as:

z ∼ N(0, Id) z1 ∼ N(0, Id1) z2 ∼ N(0, Id2)
x1 | z, z1 ∼ N(µ1 + W1z + L1z1, Ψ1)
x2 | z, z2 ∼ N(µ2 + W2z + L2z2, Ψ2) (3)

where d is the number of across-area latent variables, and d1 and d2 are the number of within-area latent variables for areas 1 and 2, respectively. The identity matrices Id ∈ Rd×d, Id1 ∈ Rd1×d1, and Id2 ∈ Rd2×d2 describe the prior covariance of each set of latent variables. Given n1 neurons from area 1, W1 ∈ Rn1×d is the loading matrix for the across-area component in area 1, and L1 ∈ Rn1×d1 is the loading matrix for the within-area component in area 1. The matrix Ψ1 ∈ Rn1×n1 is a diagonal matrix containing the independent variance of each neuron in area 1, and µ1 ∈ Rn1 represents the mean spike counts of area 1. The parameters for area 2 are defined analogously.

From Equation (3), the marginal distributions for x1 and x2 are:

x1 ∼ N(µ1, W1W1T + L1L1T + Ψ1)
x2 ∼ N(µ2, W2W2T + L2L2T + Ψ2) (4)

We can observe that pCCA-FA decomposes the covariance of neural activity in area 1 as the sum of an across-area component (W1W1T), a within-area component (L1L1T), and an independent neuron component (Ψ1). The covariance of neural activity in area 2 is decomposed into analogous across-area (W2W2T), within-area (L2L2T), and independent neuron (Ψ2) components. In the pCCA-FA model, variance shared with neurons in the other area is assigned to the across-area component (appearing in the pCCA-FA model as Wm for each area m, and corresponding to latent variable z, which involves both areas), regardless of whether that variance is also shared with neurons in the same area. Consequently, variance is assigned to the within-area component only when it is not shared with any neuron in the other area (appearing in the pCCA-FA model as Lm, and corresponding to latent variable zm, which only involves area m).

pCCA-FA can be viewed as a combination of two existing latent variable methods, namely pCCA, which finds co-fluctuation patterns that maximize correlation between two populations, and FA, which maximizes covariance between neurons within a single population. Importantly, pCCA-FA explicitly identifies both types of covariance in a single probabilistic framework. The key difference between pCCA-FA and pCCA is that pCCA-FA constrains the full covariance matrices identified by pCCA (e.g., Ψ1 in Equation (2)) as low-rank plus diagonal (e.g., L1L1T + Ψ1 in Equation (4), see Supplementary Fig. 4). It follows that pCCA-FA has fewer parameters than pCCA, and is thus more robust in the regime of low trial counts (i.e., able to recover ground truth with fewer trials; Fig. 3e-f). With sufficient trials, pCCA-FA and pCCA identify the same across-area co-fluctuation patterns.

4.12 Fitting pCCA-FA
To fit pCCA-FA to neural activity, we define a joint vector of neural activity in both areas (x˜) and a joint vector of across- and within-area latent variables (˜z):

x˜ := [x1; x2] ˜z := [z; z1; z2]

The probabilistic model of pCCA-FA from Equation (3) can be rewritten as:

˜z ∼ N(0, Id+d1+d2)
x˜ | ˜z ∼ N(µ˜ + L˜˜z, Ψ˜) (5)

where

µ˜ = [µ1; µ2] L˜ = [W1 L1 0; W2 0 L2] Ψ˜ = [Ψ1 0; 0 Ψ2]

and µ1, µ2, W1, W2, L1, L2, Ψ1, and Ψ2 are defined as in Equation (3). The formulation in Equation (5) makes clear how we can think of pCCA-FA as a constrained and structured FA model. The zeros in the loading matrix L˜ ensure that the within-area latent variables (z1 and z2) are only used to explain activity in their respective areas, whereas across-area latent variables (z) are used to explain activity in both areas.

Based on Equation (5), the marginal distribution of the activity in both areas is:

x˜ ∼ N(µ˜, L˜L˜T + Ψ˜) (6)

We fit pCCA-FA model parameters to neural activity using an exact expectation-maximization (EM) algorithm [147]. The EM algorithm for pCCA-FA resembles that used for FA. Specifically, the M-step update equations for the parameters L˜ and Ψ˜ mirror those derived in FA for the loading matrix and noise covariance matrix. Each M-step parameter update for pCCA-FA (i.e., for parameters W1, L1, etc.) effectively applies a zeroing mask to the contents of L˜, ensuring that the blocked structure of the loading matrix is maintained. Additionally, there is no constraint on the relationship between the across- and within-area co-fluctuation patterns (e.g., they are not constrained to be orthogonal).

The columns of each loading matrix (Wm and Lm for each area m) as estimated using the EM algorithm are not ordered. As a post-processing step for all analyses, we applied the singular value decomposition to transform each loading matrix, such that the first column contained the co-fluctuation pattern that explained the most shared variance, the second column explained the second-most shared variance, etc. For across-area co-fluctuation patterns, it may also be desirable to order the columns according to shared correlation (as in CCA, see below).

For each session, we jointly chose the dimensionalities for across- (d) and within-area (d1 and d2) latent variables using 10-fold cross-validation and selected the dimensionalities which maximized the cross-validated data likelihood P(x˜) (see Equation (6)). We tested integer dimensionality values between 1 and 15, inclusive, for both across- and within-area latent variables for each cross-validation run. After across- and within-area dimensionalities were selected, a final pCCA-FA model was fit with the specified dimensionalities using all trials. The parameters of this final model were used in all analyses.

4.13 Estimating latent variables in pCCA-FA
To estimate the across- and within-area latent variables in pCCA-FA (e.g., for use in predicting pupil diameter in Fig. 6), we computed the posterior mean of latent variables:

E[˜z | x˜] = L˜T (L˜L˜T + Ψ˜)−1 (x˜ − µ˜) (7)

where the first d entries of E[˜z | x˜] are E[z | x1, x2] (across-area latent variables), the next d1 entries are E[z1 | x1] (within-area latent variables for area 1), and the final d2 entries are E[z2 | x2] (within-area latent variables for area 2).

4.14 Computation of canonical variables
The pCCA-FA parameters allow for the computation of the canonical variables analogous to those found for CCA (Supplementary Fig. 6a). Following ref. [74], the canonical directions and correlations can be obtained applying the singular value decomposition on the cross-correlation matrix K ∈ Rn×n:

K = Σ11−1/2 Σ12 Σ22−1/2

where n is the total number of neurons (i.e., n1 + n2),

Σ11 = W1W1T + L1L1T + Ψ1
Σ22 = W2W2T + L2L2T + Ψ2
Σ12 = W1W2T

are the blockwise entries of L˜L˜T + Ψ˜. The canonical pairs of directions (u1i ∈ Rn1, u2i ∈ Rn2) for i = 1, 2, ..., d, are defined as:

(u1i, u2i) = (Σ11−1/2v1i, Σ22−1/2v2i),

where (v1i ∈ Rn1, v2i ∈ Rn2) are a pair of left-right singular vectors of K. The corresponding singular value ρi is the canonical correlation.

Note that the Σ matrices are defined in terms of the estimated pCCA-FA parameters rather than the sample covariance matrices (which is what CCA and pCCA use). This is the regularization that pCCA-FA provides that enables it to perform well in a limited data regime (Fig. 3e-f).

4.15 Measuring %sv and dshared
We defined two metrics to summarize the outputs of pCCA-FA: the percentage of variance shared among neurons (%sv) and shared dimensionality (dshared), which were originally introduced for FA [25, 72] and adapted here for pCCA-FA.

4.15.1 %sv
We used %sv to measure the prominence of (i.e., the proportion of variance explained by) the across- or within-area components. For a given pCCA-FA model, there are four instances of %sv: an across-area and a within-area instance, for each of the two areas. The across-area %sv quantifies the percentage of each neuron’s variance that was explained by across-area latent variables (i.e., the percentage of each neuron’s activity that was shared with at least one other neuron in the other area). More specifically:

Across-area %sv for neuron k in area m = (wmk wmkT) / (wmk wmkT + lmk lmkT + ψmk)

where wmk is the kth row of the across-area loading matrix for area m, lmk is the kth row of the within-area loading matrix for area m, and ψmk is the independent variance for the kth neuron in area m (i.e., the kth entry of the diagonal of Ψm). The numerator represents the across-area variance for neuron k in area m, and the denominator represents the total variance for that neuron. In our analyses of %sv, we averaged over neurons in a brain area to obtain a single across-area %sv metric for each area.

Similarly, the within-area %sv quantifies the percentage of each neuron’s variance that was explained by within-area latent variables (i.e., the percentage of each neuron’s activity that was shared with at least one other neuron in the same area, but was not shared with any neurons in the other area). More specifically:

Within-area %sv for neuron k in area m = (lmk lmkT) / (wmk wmkT + lmk lmkT + ψmk)

4.15.2 dshared
We used the metric, dshared, to measure the complexity of across- and within-area interactions. For a given pCCA-FA model, there are four instances of dshared, an across-area and a within-area instance, for each of the two areas. Note that dshared measures the dimensionality of interactions shared among neurons and excludes activity related to each neuron’s independent variance. We measured dshared in an analogous way to ref. [25], but adapted for across- and within-area dimensionalities.

Briefly, the first step in measuring across-area dshared is to fit the loading matrices Wm with d* across-area latent variables, where d* is the number of latent variables that maximizes the cross-validated data likelihood. Across-area dshared is defined for each area m as the minimum number of latent variables needed to explain at least 95% of the across-area shared covariance in WmWmT, where dshared ≤ d*. A 95% threshold was used in order to focus on latent variables that explain a large amount of shared variance. Note that although there is a single set of d* across-area latent variables, it is possible to obtain a different value of across-area dshared for each area.

We analogously defined within-area dshared for each area using the procedure described above by replacing WmWmT with LmLmT.

4.16 Simulated data generation
To validate the pCCA-FA model (Fig. 3e-f), as well as generate chance distributions (Fig. 5d), we generated simulated data xm for each area m from the pCCA-FA model according to Equation (4). This involved defining ground truth model parameters µm, Wm, Lm, and Ψm for each area m.

For model validation in Fig. 3e-f, we defined these parameters in the following way. We generated random loading matrices Wm and Lm for each area m using independent draws from a standard normal distribution. We scaled the %sv in each loading matrix by adjusting the singular values of the across- and within-area loading matrices (Wm and Lm) relative to each other to achieve the desired ratios (≈ 15% across-area and ≈ 10% within-area). The independent variance values (i.e., entries of Ψm) were then scaled to achieve a level of 75%. Each simulated dataset had 30 neurons per area, 5 across-area latent variables (d = 5), and 3 within-area latent variables per area (dm = 3 for each area m). For each trial count configuration, we generated 30 sets of ground truth parameters and a single simulated dataset per set of parameters.

For Fig. 5d, we defined the ground truth model parameters to be copies of those that were estimated from each of the neural recordings, with one manipulation, described below. For the chance distributions in Fig. 5d and the simulated data in Supplementary Fig. 3, we specified the angle between across- and within-area co-fluctuation patterns (θsim). To achieve this, we held the across-area loading matrix (Wm) constant and set the first column of the within-area loading matrix (Lm) to be a vector that was θsim degrees apart from the first column of Wm. The remaining columns of the within-area loading matrix were randomly generated from a standard normal distribution, and then orthogonalized relative to the first column (i.e., the first column remained fixed and the remaining columns were set to be mutually orthogonal to it).

4.17 Fitting pCCA-FA and pCCA to simulated data
To fit pCCA-FA to simulated data for the analysis in Fig. 3e-f, we used the EM algorithm described above. Across- and within-area dimensionalities were estimated in a cross-validated manner. Candidate dimensionalities for across- and within-area ranged from 1 to 10 (inclusive), for both across- and within-area latent variables. After the across- and within-area dimensionalities were selected, a final pCCA-FA model was fit with the specified dimensionalities using all simulated trials. The parameters of this final model were used to compute %sv and dshared.

To fit pCCA, we used an analogous EM algorithm, presented in ref. [74]. We fit pCCA to the same sets of simulated data as pCCA-FA. The across-area dimensionality was estimated in a cross-validated manner. Candidate dimensionalities for the across-area latent variables ranged from 1 to 10 (inclusive). After the across-area dimensionality was selected, a final pCCA model was fit with the specified dimensionality using all simulated trials. The parameters of this final model were used to compute across-area %sv and dshared.

4.18 Pupil diameter regression
To explain pupil-related variance using the latent variables estimated by pCCA-FA, the across-area latent variables E[z | x1, x2] and the within-area latent variables E[z1 | x1] and E[z2 | x2] were estimated using Equation (7). We regressed either the across-area or the within-area (right or left PFC) latent variables onto the trial-by-trial pupil diameter using a linear regression model for each set of latent variables. We then assessed r2, the proportion of variance in the trial-by-trial pupil diameter explained by the regression model.

To compute chance distributions, we repeated the procedure above, except that we performed the regression using latent variables estimated from the neural activity from a different session than the one in which we measured pupil diameter (i.e., the latent variables from session i were used to explain variance in the pupil diameter from session j, where i ≠ j). Trials were truncated in the session with more trials to ensure equal trial numbers. Chance distributions contained 240 samples for Monkeys Pe and Sa with 16 sessions, and 72 samples for Monkey Wa with 9 sessions. We report the 95% confidence intervals of this chance distribution (Fig. 6c, gray bars, termed shuffle). An alternative method to generate the chance distribution involves performing a within-session shuffle, in which the latent variables from session i are used to explain variance in the pupil diameter from session i but with permuted trial labels. This procedure yielded a similar, but slightly more lenient, chance distribution.

To control for differences in the number of across- and within-area latent variables, we repeated the entire procedure using only one predictor (i.e., the top across- or within-area latent variable) rather than the full set of predictors (i.e., d across- or dm within-area latent variables). This yielded similar results (p = 4.3e−22; pooled across animals and brain areas, n = 82 brain areas, 2 per session; one-tailed t-test on the differences in r2 for each session).

4.19 Statistical tests and reproducibility
Data preprocessing was performed using custom code in MATLAB (Mathworks, Natick, MA). Criteria for inclusion of units are described in Methods section: Preprocessing of neural data. All analyses were performed in Python (version 3.8.18). Statistical tests (t-test) were performed using the SciPy (version 1.10.0, https://scipy.org) statistics module. Sample sizes and p-values are directly reported in the corresponding figure caption. Averaged data shown in figures are averaged across sessions (n=42) or brain areas (n=84, 2 per session), noted in each caption. For data shown in histograms, means are marked with triangles. Detailed descriptions of plotted data are provided in each caption. Code to reproduce the figures and statistics is publicly available on GitHub at https://github.com/meganmcd13/figures-dual-pfc.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_2": "This paper develops a method combining canonical correlation analysis and factor analysis and applies it to simultaneously recorded spike activity from both sides of the monkey's lateral prefrontal cortex during a delayed eye movement task, and decomposes the intertrial variability of the activity into a common component between the left and right prefrontal areas and a component within each area. By showing that the common inter-areal component corresponds well to the pupil diameter reflecting the arousal level, they demonstrated that this method can explain specific phenomena effectively.\nBoth canonical correlation analysis and factor analysis are major data analysis methods, but I appreciate the fact that this method integrates them in a straightforward manner to obtain previously unobtainable results. The description is also clear and easy to understand, and the validity of the method and results is well verified. Therefore, I have no major objections to the scientific significance of this paper, and I will limit myself to making a few not-so-serious points regarding the acceptance or rejection of this paper. However, I have added some comments at the end that are not related to the acceptance or rejection of the paper. I hope they will find them useful in their future research.\nMajor\n・The limitations of this method have not been discussed. What is supposed to be assumed by this method? What kinds of things might be missed? It seems natural, in a sense, that left and right hemisphere activity would co-vary with arousal state, but why was this not seen in a simple correlation and could only be obtained with the present method? etc. These questions need to be answered.\n・Canonical correlation analysis assumes that correlations exist. Will this method ever produce a pseudo correlation? In other words, will data from different trials or from different animals, measured under the same task conditions, produce any pseudo-positive results?\n・Is it possible to interpret the significance of the dimensions obtained? If this question is wrong, there is no need to answer it.\nMinor\n・I suppose they have performed single-unit sorting. If my understanding is correct, describe it clearly.\n・If rsc is calculated from multi-unit activity, will inter-regional correlations be obtained?\n・Is it possible to interpret the significance of the dimensions obtained? If this question is wrong, there is no need to answer it.\n・Schwartz and Goldman-Rakic's anatomical study [82] has indicated that nerve connections between the left and right prefrontal cortex are “comb-like.” The Utah array may also include areas that are directly innervated from the contralateral side and areas that are not. The correlations may then differ significantly depending on the electrodes in the Utah array.\n・It would be helpful to the reader if the recovery of the %sv and dshared in Fig. A3 had a figure that clearly explains why this method requires less data than pCCA. If it is not easy to illustrate, it is not required.\nComments\nThis method is interesting as a data analysis method and has a wide range of applications. For example, it may be used to determine the extent to which fluctuations in U.S. and Japanese stock prices are due to bilateral and domestic factors, respectively. Of course, the method could be applied to LFP (local field potential) as well as to spike activity. I have no doubt about the usefulness of this method.\nHowever, I did not find in this paper much willingness or a clear perspective to delve into the nature of how the brain interacts between the cerebral hemispheres.\nIf they wanted to examine interhemispheric interactions using the ODR task, they should have included an anti-saccade task or some other experimental paradigm that required specific information to be exchanged between the two hemispheres.\nThe bin width of 1 second does not seem appropriate either. The exchange of information between hemispheres necessary to accomplish a task should occur in a much shorter time span. I know that many neurons in the prefrontal cortex have low firing rates. If so, then the instantaneous firing rate or other similar index should be used.\nWhat information is carried from the contralateral area and how is it carried? Is it carried in different rhythms depending on the situation or is there a specific time window? Does a variety of information arrive unselected or is it carefully chosen? Some diseases, such as patients with apraxia, have problems with the consistency of information processing between hemispheres. It is the mission of brain scientists to try to get to the essence of how functional differentiation and integration between hemispheres is achieved.", "reviewer_3": "In their manuscript “Interactions across hemispheres in prefrontal cortex reflect global cognitive processing”, McDonnell et al. propose a method termed pCCA-FA (probabilistic canonical correlation analysis in combination with factor analysis) to study interactions within and between the left and right dorsolateral prefrontal cortex (also termed regions or areas here). For each neuron, pCCA-FA estimates the proportion of trial-by-trial spike count variance attributable to self, within-region, and across-region sources. The authors also apply SVD to extract the main latent variables underlying these sources. They show that this approach offers advantages over previous methods (rSC, pCCA, or FA). They applied this method to study the populations of each prefrontal cortex in rhesus monkeys trained on a memory-guided saccade task. They find that the primary latent variable representing the across-area co-fluctuation pattern is a more effective predictor of pupil diameter than the corresponding within-area latent variable. Since pupil diameter is a marker of arousal, a brain-wide process, this finding supports that pCCA-FA can find truthful and insightful patterns and sources of variability in the brain.\nI thought the approach the authors propose to study neuronal interactions is elegant and simple. pCCA-FA combines two well-established methods and, as far as I understand, is solidly constructed. They also provide an adequate justification for the need for a novel analysis tool. However, aspects related to the application of the method to the dataset and the interpretation of the results could be further refined.\nGeneral points:\n1. There is a discrepancy between the title and narrative of the manuscript. From the title, one would think that the main focus of the paper is the study of “global cognitive processing,” which is indeed addressed in Results section 2.5, where pupil diameter is analysed. However, the discussion recurrs only briefly to this topic. I would like to see more efforts to relate the findings to cognition or, alternatively, a change in the title.\n2. It is unclear to me why both the “mean pupil diameter” and the “pupillary evoked response magnitude” are analysed. The differences between these two measurements, as well as their cognitive interpretation, are not discussed.\n3. In the discussion, it is mentioned that anatomical constraints on network circuitry may affect the correlation between neurons [line 363]. However, the cited papers are about visual, not prefrontal cortex. I would sugges to include more appropriate references, e.g. Kritzer & Goldman-Rakic, 1995, Journal of Comparative Neurology; Pucak, Levitt, Lund, & Lewis, 1996, Journal of Comparative Neurology\n4. The results for the two different eye-pupil measures appear to differ slightly, particularly for monkey “Wa.” Moreover, monkey “Wa” shows overall different values for rsc, %sv, and dshared. Is there any interpretation of these differences?\nParticular points:\n[103-105] The phrasing in this sentence could be clearer. Since coordination between neurons is defined as synchronized patterns of activity, describing it as a process that causes or “results” in such patterns may appear tautological. This could also unintentionally suggest that neurons coordinate through signals other than their activity. Rewording the sentence could help avoid this ambiguity.\n[140] What statistical test was performed here? If I understood correctly, none of the tests shown in Fig. 2 directly compare the distributions of rSC within versus across regions.\nAdditionally, for monkey “Wa,” the distributions of average rSC within and across regions across all sessions do not appear very different to the naked eye.\n[248-251] Is this discrepancy the only reason that rSC is smaller than %sv? Could it be resolved simply by taking the absolute value or the squared value (as is done with %sv)?\n[316-317] The “pupillary evoked response magnitude” is measured “just after visual stimulus presentation,” which occurs before the delay period. From the current wording (“across-area latent variables were able to explain multiple aspects of pupil-linked arousal”), one might infer a causal role from the latent variable to the pupillary evoked response. However, the neuronal activity is recorded “1 second of the delay period immediately preceding the earliest possible saccade initiation cue in the session.” Therefore, it cannot be concluded that the “latent variables identified by pCCA-FA explained changes in pupil diameter,” since the change in pupil diameter precedes the neuronal activity, implying that the causal direction would be reversed.\n[331-332] I would argue that this is not, in general, “the shared variance among PFC neurons.” Here, task conditions (right or left stimulus cue, or right or left response) are pooled together. Therefore, it specifically measures the task-condition-independent variance. This distinction is significant because the PFC is primarily involved in processes that are task-dependent, where the task condition plays an important role. Thus, the analysis is capturing the shared variance of a specific sub-process within the PFC.\n[481-485] Why was the entire electrode signal removed instead of only removing the repeated units? Could this approach be removing nearby neurons that might be more connected, thereby further reducing the %sv of inter-region units?\n[500-502] If I understand correctly, the intention of this procedure is to remove condition-specific rate properties. However, this method only removes cue-specific rate variations, omitting response-specific variations. In other words, when pooling both correct and incorrect trials, for the same cue condition, there can be two overlapping modes of firing rate depending on whether the response was right or left. The paper does not specify which trials were selected (all or only correct) or what the hit rate was.\nFurthermore, it would strengthen the manuscript to demonstrate that condition-mean subtraction effectively removes condition-specific activity patterns. This could be done, for example, by training a model to classify trials by condition before and after rate manipulation, showing that condition-mean subtraction renders the conditions indistinguishable.\n[503] It is mentioned that z-scoring “yielded qualitatively similar rSC results,” but nothing is said about its effect on the pCCA-FA results. Should we expect that z-scoring would also lead to similar outcomes in pCCA-FA as in rSC?\n[818-827] Why was the null distribution of the regressions created by shuffling neuronal activity from neurons in other sessions? Why not simply shuffle the pupil diameter values within a session while using the firing rates of the same neuron but from other trials?", "reviewer_4": "This is an interesting and important paper of potentially high impact describing a new statistical technique developed by the authors that can parse components of variation in population activity patterns attributable to local interactions between neurons and inputs from neurons in another brain area.\nThis is a key challenge in deciphering the computational organization of the cerebral cortex, and it is the advantage of this technique over alternative methods. The technique could significantly advance our understanding of the interplay between local processing in cortical areas and large scale communication in cortical networks.\nNearly all of the questions I initially had about the validity of the technique were answered by results provided in Extended data. Generally, I recommend promoting as many of the highlights of the Extended results into additional panels of the main figures as possible, so that the majority of the readership sees the strong evidence validating the new technique.\nIn addition to describing and validating a new analytical technique, the paper provides important new information about the coordination of neural activity in prefrontal cortex, including the findings that correlations between brain areas (left and right hemispheres) are stronger than correlations within areas, and that the strength of within and between area neural correlations track several behavioral variables.\nMajor\n1. The paper could be strengthened by enhancing the presentation of the data validating the pCCA-FC technique. Much of the crucial information is provided in the Extended Data. This includes: (1) comparison of detected correlations in the original data with trial shuffled data (Fig. A1), (2) relation of noise correlation to signal correlation (Fig. A2), (3) recovery of ground truth parameters from simulated data (Fig. A3), and (4) evidence that the pCCA-FC technique outperforms alternative methods (FIg. A3). The paper could be strengthened if some of the key panels from these figures were promoted to the main figures of the paper (for example adding them to Fig. 3). This would strengthen the case that the analysis developed by the authors measures real neural correlations, successfully parses variation within and between cortical areas, and outperforms existing methods\n2. Related to the above, trial shuffling is utilized in the extended data to validate removal of slow sources of correlation over trials (Fig. A1). The current data could be strengthened by showing that trial shuffling destroyed detected neural correlations at faster time scales after removal of the slower components of correlation. This would substantiate that detected neural correlations depended on the simultaneity of the neural signals and were not artifacts of the analysis. The amount of correlation detected in the trial shuffled data (in some cases ~10%, Fig. 1Ad) was not trivial in relation to the magnitude of the correlation in the original data (~10-15%, Fig. 4B). It would therefore be important to show that the neural correlations detected by this method significantly exceed the level of correlation seen after trial-shuffling.\n3. A key claim/finding is that the new analysis can parse correlations driven by local interactions between neurons from correlations driven by input from other brain areas. If this is true, differences in %sv attributable to interactions within and between cortical areas should be greater in the original data than after shuffling neurons between brain areas (randomly changing right/left hemisphere labels). Demonstrating that with an additional bootstrap test if feasible would strengthen the main claim of the paper.\n4. More %sv was attributed to interactions between than within brain areas. The models derived to estimate %sv between areas had more free parameters (a greater number of latent parameters). It will be important for the main findings of the paper to argue/show that the difference in %sv attributed to local and long tract interactions did not reflect differences in the number of latent parameters in the models used to fit the two types of interactions.\n5. The analysis is built to detect noise correlations, e.g. correlations in neural signals after modulation of activity due to task performance has been subtracted out. Relations to behavior (pupil diameter, spatial tuning) are therefore indirect. Can the technique be adapted to detect correlations in neural signal (task-related modulations in firing rate)? This would enhance links between detected correlations and behavior. Some description of how this might be done in future studies in the Discussion could be useful to understand the potential applicability of the technique to different questions.\n6. It would be useful to clarify the link between detected neural correlations and the physiological processes that are likely to generate them, namely synaptic interactions between neurons. In particular, how do the authors envision that synaptic interactions between neurons could drive correlations in firing rate at the trial level that are captured by their analysis? For example, correlations in delay period firing rate over trials could reflect net excitatory and inhibitory synaptic interactions between neurons integrated over the 1 sec delay period, etc.\n7. Related to point 6 above, the authors hypothesize that inter-area correlation is driven by a global neuromodulatory signal (such as NE) potentially linked to arousal. If the signal is global, why would it preferentially strengthen correlations between neurons in different brain areas in comparison to correlations between neurons in the same area? The logic here needs additional clarification.\n8. Also related to point 6 above, it is thought that more synapses onto cortical pyramidal neurons derive from neighboring neurons in the same brain area than from neurons in distant areas. If synapses from recurrent local circuits do in fact provide the bulk of synaptic input, why then do inputs from remote brain areas drive stronger correlations? It could be useful to consider the correlation results in the context of what is known about the architecture of local versus long-tract connections, and to reconcile apparent discrepancies between the relative strength of connections and the relative strength of correlations.\nMinor\n1. The Methods describe training a neural network to distinguish neural action potential waveforms from noise. Were action potential waveforms from different neurons discriminated? If not, it would be important to state that the analysis was conducted on multi-neuron firing rates.\n2. The paper does an excellent job of introducing the statistical concepts behind the technique. Any additional walk-throughs of the main concepts that the authors can provide would be to the good, to clarify the relationship between the loading matrices, latent variables, and population activity vectors, and how they combine to generate the co-fluctuation in firing patterns illustrated schematically in Fig. 3b.", "reviewer_1": "Response to reviewers\nInteractions across hemispheres in prefrontal cortex reflect global\ncognitive processing\nMegan E. McDonnell*, Akash Umakantha*, Ryan C. Williamson*,\nMatthew A. Smith†, & Byron M. Yu†\n* and † indicate equal contribution\nWe would like to thank the reviewers for their thoughtful and constructive comments. We have\nperformed new control analyses and made extensive updates to the figures, supplements, and\ntext throughout the manuscript, especially in the Discussion. We believe the manuscript is now\nstronger and clearer as a result.\nThe main revisions are as follows:\n● We have significantly revised the Discussion as suggested by multiple reviewers.\nSpecifically, we have incorporated feedback to address limitations and applications of\npCCA-FA (Reviewers 2 & 3), and included more discussion of how the brain may interact\nbetween hemispheres at the cellular and circuit levels (Reviewers 2 & 4).\n● We have promoted results from Supplementary Figs. 2 and 3 (numbers reflect figure\norder at initial submission) into the main figures to better highlight the validation of\npCCA-FA and the link between r and signal tuning, as suggested by Reviewer 4.\nsc\n○ Supplementary Fig. 2 (previous number) has been promoted to Fig. 2e.\n○ Supplementary Fig. 3c-d (previous number) have been promoted to Fig. 3e-f.\n● We have performed additional analyses requested by the reviewers and added 3 new\nsupplementary figures (now Supplementary Figs. 5, 7, and 8).\n○ Supplementary Fig. 5 shows a shuffle control in which randomizing the\nhemisphere labels destroys within-area %sv in the PFC recordings, as suggested\nby Reviewer 4.\n○ Supplementary Fig. 7 shows a simulation that considers the effect of\nsubsampling different subsets of neurons from session to session on the\nvariation of θ_estimated that we report in Fig. 5c, as suggested by Reviewer 1.\n○ Supplementary Fig. 8 is devoted to additional analyses of the pupil. We illustrate\nour preprocessing steps, definition of event-related (also known as evoked) pupil\ndiameter, and regression onto event-related pupil diameter (previously Fig. 6d).\nThis reorganization was in response to multiple reviewers (1, 3 & 4).\nAccompanying this change were updates to the text surrounding the main text\npupil figure (Fig. 6).\n● Beyond the changes listed above, we have added 1 new main figure panel, 2 new\nsupplementary figure panels, and made updates to 1 supplementary figure panel.\n○ Fig. 3b: Reviewer 4 suggested that we clarify the relationship between the\nloading matrices, latent variables, and population activity vectors. To address\nthis, we added a schematic illustrating the use of co-fluctuation patterns to map\nlatent variables to neural activity.\n○ Supplementary Fig. 1e: Reviewer 4 suggested that we include an additional\nshuffle control to Supplementary Fig. 1 which expands on our result that shuffling\nthe trial labels destroys detected neuronal correlations. To address this, we have\nadded a panel analogous to Supplementary Fig. 1f (previously 1e) which uses\nthe shuffle control to demonstrate that the level of correlation in the PFC\nrecordings is higher than chance.\n○ Supplementary Fig. 3a: Reviewer 1 suggested that we adjust how we illustrate\nthe accuracy of pCCA-FA when the angle between across- and within-area\nco-fluctuation patterns is small. To address this, we have remade Supplementary\nFig. 3a to more clearly show the scale of estimation errors made by pCCA-FA,\nas well as adjusted the accompanying text to make more precise claims.\n○ Supplementary Fig. 4b: Reviewer 2 suggested that we illustrate the number of\nparameters in pCCA-FA relative to pCCA. In response to this we have added a\nnew panel to Supplementary Fig. 4 which graphically illustrates the parameter\ncounts using simulated data.\nBelow, we have provided a point-by-point response (in pink) to each reviewer comment (in\nblack). All line numbers and figure numbers refer to the updated manuscript unless specified.\nReviewer #1 (Remarks to the Author):\nThe manuscript “Interactions across hemispheres in prefrontal cortex reflect global cognitive\nprocessing” by McDonnell, Umakantha, Williamson and colleagues shows that a large\nproportion of neural variability is shared between left and right prefrontal regions. They\nfurther show that this shared variability is linked to an arousal process, expressed in pupil\ndilation. To reach these conclusions they analyzed data using a novel method, pCCA-FA,\nwhich allows the partition of across- and within-area variance.\nOverall, the manuscript is well-written and clear. The methods seem appropriate and well-\njustified. pCCA-FA will become an important tool for future research. The research findings are\nnovel and interesting. Below I list some clarifications that I think will make the\nmanuscript more clear to the reader.\nWe thank the reviewer for their positive comments.\nMajor:\nL518: I have a doubt about the method used to separate slow and fast components. In line 518\nthe authors write:\n“We computed the slow component for each session using a centered boxcar filter of length 25\ntrials"}, "author_response": "We would like to thank the reviewers for their thoughtful and constructive comments. We have performed new control analyses and made extensive updates to the figures, supplements, and text throughout the manuscript, especially in the Discussion. We believe the manuscript is now stronger and clearer as a result.\nThe main revisions are as follows:\n● We have significantly revised the Discussion as suggested by multiple reviewers.\nSpecifically, we have incorporated feedback to address limitations and applications of pCCA-FA (Reviewers 2 & 3), and included more discussion of how the brain may interact between hemispheres at the cellular and circuit levels (Reviewers 2 & 4).\n● We have promoted results from Supplementary Figs. 2 and 3 (numbers reflect figure order at initial submission) into the main figures to better highlight the validation of pCCA-FA and the link between r and signal tuning, as suggested by Reviewer 4.\nsc\n○ Supplementary Fig. 2 (previous number) has been promoted to Fig. 2e.\n○ Supplementary Fig. 3c-d (previous number) have been promoted to Fig. 3e-f.\n● We have performed additional analyses requested by the reviewers and added 3 new supplementary figures (now Supplementary Figs. 5, 7, and 8).\n○ Supplementary Fig. 5 shows a shuffle control in which randomizing the hemisphere labels destroys within-area %sv in the PFC recordings, as suggested by Reviewer 4.\n○ Supplementary Fig. 7 shows a simulation that considers the effect of subsampling different subsets of neurons from session to session on the variation of θ_estimated that we report in Fig. 5c, as suggested by Reviewer 1.\n○ Supplementary Fig. 8 is devoted to additional analyses of the pupil. We illustrate our preprocessing steps, definition of event-related (also known as evoked) pupil diameter, and regression onto event-related pupil diameter (previously Fig. 6d). This reorganization was in response to multiple reviewers (1, 3 & 4). Accompanying this change were updates to the text surrounding the main text pupil figure (Fig. 6).\n● Beyond the changes listed above, we have added 1 new main figure panel, 2 new supplementary figure panels, and made updates to 1 supplementary figure panel.\n○ Fig. 3b: Reviewer 4 suggested that we clarify the relationship between the loading matrices, latent variables, and population activity vectors. To address this, we added a schematic illustrating the use of co-fluctuation patterns to map latent variables to neural activity.\n○ Supplementary Fig. 1e: Reviewer 4 suggested that we include an additional shuffle control to Supplementary Fig. 1 which expands on our result that shuffling the trial labels destroys detected neuronal correlations. To address this, we have added a panel analogous to Supplementary Fig. 1f (previously 1e) which uses the shuffle control to demonstrate that the level of correlation in the PFC recordings is higher than chance.\n○ Supplementary Fig. 3a: Reviewer 1 suggested that we adjust how we illustrate the accuracy of pCCA-FA when the angle between across- and within-area co-fluctuation patterns is small. To address this, we have remade Supplementary Fig. 3a to more clearly show the scale of estimation errors made by pCCA-FA, as well as adjusted the accompanying text to make more precise claims.\n○ Supplementary Fig. 4b: Reviewer 2 suggested that we illustrate the number of parameters in pCCA-FA relative to pCCA. In response to this we have added a new panel to Supplementary Fig. 4 which graphically illustrates the parameter counts using simulated data.\nBelow, we have provided a point-by-point response (in pink) to each reviewer comment (in black). All line numbers and figure numbers refer to the updated manuscript unless specified.\nReviewer #1 (Remarks to the Author):\nThe manuscript “Interactions across hemispheres in prefrontal cortex reflect global cognitive processing” by McDonnell, Umakantha, Williamson and colleagues shows that a large proportion of neural variability is shared between left and right prefrontal regions. They further show that this shared variability is linked to an arousal process, expressed in pupil dilation. To reach these conclusions they analyzed data using a novel method, pCCA-FA, which allows the partition of across- and within-area variance.\nOverall, the manuscript is well-written and clear. The methods seem appropriate and well-justified. pCCA-FA will become an important tool for future research. The research findings are novel and interesting. Below I list some clarifications that I think will make the manuscript more clear to the reader.\nWe thank the reviewer for their positive comments.\nMajor:\nL518: I have a doubt about the method used to separate slow and fast components. In line 518 the authors write:\n“We computed the slow component for each session using a centered boxcar filter of length 25 trials, computed after removing target information as described above. We then subtracted this component from the residual spike counts to remove slow-timescale trends that could have induced spurious correlations”\nCan the authors provide some validation to the claim that this removes the slow component? I ask because I have seen recordings where non-biological technical issues (such as clustering one neuron into two different spikes due to slow drift, or simply electrode movements that make you lose a spike half-way through the recording session) may lead to any number of non-stationarities in the data that would not be removable with the 25-trial filter that the authors employ. As an example, below I show a reanalysis of Panichello and Moore Nature (2024), who used neuropixels from the same region (area 8) using the same task (delay saccade task) as in the manuscript. Here I am plotting the average delay period activity of several neurons (each plot is one neuron) from one recording session (each point is the mean delay activity of a trial, and the blue ticks can be ignored). The x-axis is the trial number. As can be seen, if the method employed in this manuscript had been applied to this dataset, the 25-trial filter would not have been effective for many of the cells.\nWe appreciate the reviewer’s concern and agree that a boxcar filter is not suitable to remove many types of artifactual non-stationarities in neural recordings. Before applying the boxcar filter, we performed a preprocessing step to remove electrodes on the Utah array that exhibited activity non-stationarities that were likely non-biological (e.g., cells appearing or disappearing during a recording session), as the reviewer illustrates. We have clarified this preprocessing step in the Methods section. In short, we removed electrodes that had a high Fano Factor (>10) from further analysis. This successfully excluded electrodes similar to those the reviewer plotted above. As an example, we provide an analogous plot of example electrodes from Monkey Sa that were removed by this cutoff, and example electrodes of those that remained as part of our analyses. On the left, we show a histogram of the Fano Factor of all electrodes in all sessions recorded for Monkey Sa. The FF=10 cutoff is marked with a black dashed line. On the right, we provide 3 example electrodes from one session that were removed for having a high Fano Factor (each marked on the left with a red dashed line), and 5 example electrodes that were kept (each marked on the left with a green dashed line).\nThe boxcar filter was then applied to electrodes that remained after the preprocessing (e.g., electrodes indicated with the green arrow in the plot above), and this served to remove what we believe to be true (i.e., biological) neural slow drift (as described in Cowley, et al. 2020). We verify that using this method destroys spurious correlations by performing the flip control in Supplementary Fig. 1. We updated the text in the Methods section to clarify that boxcar filtering is only applied to electrodes kept after the removal of artifactual electrodes.\nChanges:\n● L546 in Methods section “Preprocessing of neural activity”:\n○ Addition: “This served to eliminate electrodes contaminated by recording artifacts, or units that were not present for the entirety of a recording session.”\n● L611 in Methods section “Separation of slow and fast components”:\n○ Original: “We computed the slow component for each session using a centered boxcar filter of length 25 trials, computed after removing target information as described above.”\n○ Updated: “We computed the slow component for each session using a centered boxcar filter of length 25 trials, computed for each remaining electrode after the preprocessing described above.”\nL220: “We verified that pCCA-FA accurately identified the ground truth %sv and d_shared using simulated data (Extended Data Fig. A3a-b and Extended Data Fig. A4)”. Extended Data Fig. A4 shows that the recovery of %sv depends on angle between across- and within-area co-fluctuation patterns. How does this show that pCCA-FA accurately identifies the ground truth? If the angle is low in the neural data then pCCA-FA would not be very accurate right?\nWe agree that the statement on L220 (original number) was too specific in its claim. We meant to say that pCCA-FA accurately identified the ground truth %sv and d_shared in simulated data using realistic numbers of trials and parameters. This is supported by the following two findings. First, in Supplementary Fig. 3 (previously Extended Data Fig. A4), even when the angle between the across- and within-area co-fluctuation patterns (i.e., θ_sim) was low, the error in recovered values of %sv was within 5%sv of the ground truth values, which were ~15% across-area %sv and ~10% within-area %sv. Second, although across- and within-area %sv were not always accurately estimated when the angle (θ_sim) is low, the total %sv (across + within) is well-estimated. To incorporate these points, we have adjusted the layout of Supplementary Fig. 3a (previously Extended Data Fig. A4a) to emphasize that the estimated values were close to the ground truth even when the angle (θ_sim) was low. Further, we updated the text referencing this figure to make a less strict claim.\nChanges:\n● L232 in Results for Fig. 3:\n○ Original: “We verified that pCCA-FA accurately identified the ground truth %sv and d_shared using simulated data.”\n○ Updated: “We verified that pCCA-FA accurately identified the ground truth %sv and d_shared using simulated data with properties chosen to be similar to our neural recordings (Supplementary Figs. 2-3).”\n● Supplementary Fig. 3a (previously Extended Data Fig. A4a):\n○ Instead of plotting estimation error in within-area %sv on the horizontal axis against estimation error in across-area %sv on the vertical axis, we now plot the ground truth %sv on the horizontal axis against the estimated %sv on the vertical axis. This layout better demonstrates that the estimated values of %sv for both across- (pink) and within-area (blue) %sv lie close to the ground truth values of %sv (points lie near diagonal). We also added the black point indicating total shared variance (across-area plus within-area %sv), which was well-estimated (i.e., close to the diagonal) in all cases.\n● Caption:\n○ Addition: “Although the individual across- (pink) and within-area (light and dark blue for areas 1 and 2, respectively) points sometimes deviated from the ground truth (indicated by their vertical offset from the dashed line), the total %sv (i.e., across-area plus within-area %sv, indicated by black point) was well estimated (close to dashed line). Thus, the estimation errors of across- and within-area %sv involved overestimates of one type of %sv (across- or within-area) with corresponding underestimates of the other, thereby providing an accurate estimate of the total ground truth level of %sv.”\nL1081 (Figure 5): “For this example session, there was a high diversity of weights, inconsistent with a global change in responsivity.” I don’t understand why this is inconsistent. Couldn’t different neurons respond differently to gain modulation? For example, if inhibitory interneurons increase their gain, pyramidal cells may decrease (or vice versa). Or is the assumption that a global change must occur in all/most cells? I think that it would be good to provide references to prior observations of the types of global changes that are being assessed here, assuming that this has been observed in the literature. If this hasn’t been observed, I suggest removing this interpretation since it would be a strawman argument.\nWe thank the reviewer for this pointing this out. In hindsight, we realized that the use of the word “responsivity” could be confusing, as the reviewer indicates. To address this, we have removed the term from the location indicated (Fig. 5) as well as from the main text. In its place, we have elaborated that we meant to assess the degree to which the neurons shared a similar signed weight onto the latent variables.\nChanges:\n● L1315 in Fig. 5 caption:\n○ Original: “For this example session, there was a high diversity of weights, inconsistent with a global change in responsivity.”\n○ Updated: “For this example session, there was a high diversity of weights, many with a positive sign and many with a negative sign.”\n● L282 in Results for Fig. 5:\n○ Original: “We first asked whether the across- or within-area co-fluctuation patterns involved the entire populations increasing or decreasing their activity together, which might be the case for a global change in responsivity.”\n○ Updated: “There are many possible population-wide co-fluctuations. Perhaps the simplest example is all neurons sharing the same-signed loading onto a latent variable, such that the population would all increase their activity together or decrease their activity together. We examined whether this was the case in our PFC recordings.”\nL290: “Together, these results indicated that across- and within-area components involved different co-fluctuation patterns that varied in their similarity from session to session.” I am not convinced that this is the best interpretation of the results. To me, it looks like you have a true angle that is close to 90, but the noisiness of the estimate across sessions leads to lower estimated angles. If this is true, the sessions with the most numbers of cells/trials should be closer to 90 degrees, and sessions with fewer cells/trials should have angles below 60. Could you take a look at this? It’s just hard to think that the same animals are employing different mechanisms in different sessions.\nWe appreciate the reviewer bringing this up, as we wondered about this as well. In fact, while preparing the original submission, we spent several months performing different analyses to understand the variation in θ_estimated across sessions, which ultimately led to Fig. 5d. To summarize this effort, we considered two categories of explanation.\nFirst, as the reviewer suggests, we determined that neither the number of trials nor the number of neurons in the session had an effect on the estimated angle (θ_estimated). To show this, we plotted the value of θ_estimated on the vertical axis against the number of trials (top), number of neurons in the left hemisphere (middle), and number of neurons in the right hemisphere (bottom). In all cases, there are points in the top left of the scatter plot, indicating that a low number of neurons or trials yielded a high θ_estimated (e.g., circled points). There did not appear to be a consistent relationship between θ_estimated and the number of trials or number of neurons. To incorporate this result into the analysis in Fig. 5d, we designed the chance distributions such that each recording session had a matched sample in the distributions obtained for the chance distributions (i.e., matched neuron count, trial count, etc.). We agree that this was not made clear in the main text, so we added a statement on the matched parameter count in the description of Fig. 5d.\nSecond, we considered the related question of whether a fixed angle between the across- and within-area components may exist in PFC as a whole, but our estimation of this fixed angle was noisy due to sampling a different subset of neurons on each recording session. To determine whether this was a plausible explanation for the width of the distribution of θ_estimated we observed in Fig. 5c, we performed a new simulation for which we have added a new supplementary figure (Supplementary Fig. 7). We found that subsampling from a larger fixed population of neurons (with angle θ_full) yielded some variation in the resulting angles of the subsampled population (top). However, in comparing the recovered values of θ_estimated (bottom) to those in Fig. 5d, we concluded that this variation was insufficient to explain the width of the distribution of θ_estimated in Fig. 5d. Instead, the results in Fig. 5 and Supplementary Fig. 7 imply the existence of some inherent variability in the angle between the across- and within-area components from session to session.\nChanges:\n● L303 in Results for Fig. 5:\n○ Original: “We simulated activity in each area with a specified value of θ that sim was otherwise matched to the neural recordings (see Methods).”\n○ Updated: “We simulated activity in each area with a specified value of θ that sim was otherwise matched to the neural recordings in terms of the number of neurons, latent variables, and trials (see Methods).”\n● L311 in Results for Fig. 5:\n○ Addition: “Furthermore, we determined that the observed variation in θ_estimated did not arise solely from sampling different subsets of neurons from session to session, although it still likely contributed (Supplementary Fig. 7).”\n● Supplementary Fig. 7:\n● Caption:\n“In Fig. 5c, we observed a distribution of angles between the top across- and within-area co-fluctuation patterns (θ_estimated). We considered two possible explanations for this variation across sessions. First, in Fig. 5d, we considered whether the variation of θ_estimated was due to fitting pCCA-FA using a finite number of noisy simulated data samples. Second, we consider here the possibility that the variation arose from sampling different subsets of neurons on each recording session from a larger fixed population.\nTo simulate the second scenario, we defined a large number of neurons (10,000) representing the entirety of each simulated brain area (e.g., PFC). For the 10,000 simulated neurons in each brain area, we generated random across- and within-area loading matrices W_m and L_m for each area m using independent draws from a standard normal distribution, and scaled their columns (d=5 across-area latent variables, and d_m=3 within-area latent variables for each area m) to ensure that the ground truth value of across-area %sv was ~15% and within-area %sv was ~10%, as in Supplementary Fig. 2. We then fixed the angle between the top across- and within-area co-fluctuation pattern, θ_full, to be either 60° or 90° (see Methods).\nRecordings with a Utah array cover only a fraction of the brain area and yield a partially overlapping subset of neurons on each recording session. To represent the pool of neurons from which we sample on each recording session, we randomly selected a subset of neurons (n_sub1=100; n_sub2=100) in each simulated brain area. Then, for each of 100 independent simulations, we randomly subsampled n_1=30 and n_2=30 neurons from the fixed pool of n_sub1 and n_sub2 neurons as the neurons “recorded” in that simulated session. To quantify the variation in angle that arises from sampling a different subset of neurons on each recording session, for each simulation, we computed the angle between the top across- and within-area co-fluctuation patterns (θ_sim) using only the n_1 and n_2 “recorded” neurons. Looking across the 100 simulations yields a distribution of θ_sim for each value of θ_full (top panel, mean of each distribution marked with a triangle). For visual clarity, we plotted the θ_sim values from only one simulated brain area (both areas yielded similar results). The width of each of these distributions (100 values in each distribution) represents the variation in angle that arises from sampling a different subset of neurons on each recording session.\nFor each simulation (corresponding to one value of θ_sim above), we generated 1,000 trials of simulated data and fit a pCCA-FA model using the ground truth values of across- and within-area dimensionality, as in Fig. 5d. Using the across- and within-area loading matrices estimated by pCCA-FA, we computed one value of θ_estimated per simulated brain area. Looking across the 100 simulations yields a distribution of θ_estimated for each value of θ_full (bottom panel, mean of each distribution marked with a triangle). As in the top panel, we plotted the θ_estimated values from one simulated brain area (both areas yielded similar results). Note that each element of the histogram of θ_estimated here corresponds to one element of the histogram of θ_sim in the top panel. The width of each of these distributions (100 values in each distribution) represents the variation in estimated angle arising both from sampling a different subset of neurons on each recording session and sampling a finite number of noisy data samples.\nInterestingly, the distributions of θ_estimated in the lower panel are wider than those reported in Fig. 5d for the same values of θ_sim (for θ_full = 60, p<0.0001, for θ_full = 90, p<0.0001; Levene's tests comparing equality of variance for the distribution of θ_estimated in the lower panel to the corresponding distribution in Fig. 5d). This indicates that some of the breadth of the distribution of θ_estimated from our neural recordings in Fig. 5c may be attributable to sampling a different subset of neurons on each recording session. These results, together with those in Fig. 5d, indicate that sampling a different subset of neurons on each recording session and fitting pCCA-FA to a finite number of noisy samples cannot fully account for the width of the distribution of θ_estimated that we observed in the neural recordings. This supports the conclusion that the across- and within-area components involved co-fluctuation patterns that varied in similarity from session to session.”\nL320: Are the results in Figure 6d independent of those in 6c? Or could 6d be explained by 6c? I was not expecting the global signal to correlate with stimulus-evoked pupillary constriction (which is a reflex mediated by the midbrain), so I am wondering whether alertness signals may interact with or fully explain the results shown in Figure 6d.\nWe thank the reviewer for this point, and agree that the link between pupillary evoked response in Fig. 6d (previous number) and the trial-by-trial pupil diameter in Fig. 6c was not fully explained. In this work, the event-related pupil diameter (or evoked response) and trial-by-trial pupil diameter are measured during different windows of each trial, and have a different baseline subtraction. The trial-by-trial pupil diameter is measured during the delay period of each trial, while the event-related pupil diameter is measured just after target presentation (manually selected on each session, as shown in Supplementary Fig. 8f-h). Thus, Fig. 6d (previous number) is not trivially explained by Fig. 6c. While the pupillary evoked response is indeed reflexive, there are studies that show cognitive modulation of evoked response magnitude (reviewed by Mathot & Van der Stigchel 2015 and Ebitz & Moore 2019). We believe such cognitive modulation is why the across-area component is related to the pupillary evoked response. In order to better motivate the concepts of event-related pupil diameter, we moved Fig. 6d to a new supplementary figure devoted to the pupil (Supplementary Fig. 8). There, we present additional detail on how the pupil diameter is preprocessed, how event-related pupil diameter is defined, and its relationship to the trial-by-trial pupil diameter (Fig. 6c). In reworking the event-related pupil diameter analyses, we also adjusted the preprocessing of the pupil diameter. We adjusted the window in which we measured the trial-by-trial pupil diameter to further differentiate it from the event-related pupil diameter. The trial-by-trial pupil diameter is now exclusively measured during the delay period of each trial, from target offset to saccade initiation cue. We also added a filter for pupil data quality, which excluded one session from our analyses of the pupil (Fig. 6 and Supplementary Fig. 8). We have updated the text and Methods to reflect these changes and better motivate the measurement of event-related pupil diameter.\nChanges:\n● L321 in Results for Fig. 6:\n○ Addition: “In addition, arousal state has been linked to changes in both the tonic diameter of the pupil [79,80] and event-related changes in pupil diameter [81,82].”\n● L343 in Results for Fig. 6:\n○ Original: “Next, we asked whether the across- and within-area latent variables could explain variance in the pupillary evoked response magnitude, defined as the fast-timescale change in pupil diameter just after visual stimulus presentation (see Methods). We found that the across-area latent variables also explained more variance in the pupillary evoked response than the within-area latent variables (right PFC: p= 0.0001, left PFC: p<0.0001; Fig. 6d).”\n○ Updated: “Furthermore, we found that the across-area latent variables explained variance in the event-related changes in pupil diameter observed just after target onset (Supplementary Fig. 8; p = 0.0001, pooled across animals and brain areas; paired one-tailed t-tests).”\n● L621 in Methods section “Preprocessing of pupil diameter”:\n○ Original: “We measured pupil diameter in the time window starting 200ms before target onset until the saccade initiation cue on every trial.”\n○ Updated: “On each trial, we measured the pupil diameter during the delay period, from target offset until the saccade initiation cue.”\n● L632 in Methods section “Preprocessing of pupil diameter”:\n○ Addition: “To ensure pupil data quality, we computed the variance across trials of the trial-by-trial pupil diameter. We then removed outlier sessions that exceeded a variance threshold of two times the interquartile range of the variance values across sessions. This was performed for each animal individually. This procedure excluded one session from Monkey Wa from analyses of the pupil (Fig. 6, Supplementary Fig. 8), resulting in n=9 sessions for analysis.”\n● Supplementary Fig. 8:\n● Caption:\n“In Fig. 6, we investigated whether the across- or within-area components of neural activity were related to the trial-by-trial pupil diameter measured during the delay period (relative to the ongoing 25-trial average, see Methods). This assessed trial-to-trial fluctuations in pupil diameter relative to the ongoing mean (with a baseline subtraction that was matched to our neural data processing), and did not explicitly assess event-related changes in the pupil locked to saccade target onset. Measures of tonic pupil size have been linked to arousal and other cognitive states [10,11]. However, cognitive modulation has also been observed in evoked or event-related pupillary dynamics [12,13]. Here, we investigated whether the across- and within-area components of neural activity were related to event-related pupil diameter and compared that to the delay period pupil measurements presented in Fig. 6.\na, While animals performed the task, we monitored the raw pupil size during each trial's delay period over the course of the session (light gray). The size of the pupil slowly fluctuated over the course of the session [1], in this example session drifting upwards. As with our preprocessing of neural activity, we estimated the slowly varying component of pupil diameter by computing a moving average of 25 trials (dark gray, see Methods).\nb, We removed the estimated slow component of pupil diameter (top) from the raw pupil diameter. This yielded the residual “trial-by-trial” pupil diameter, one value per trial (bottom). In Fig. 6, we found that the trial-by-trial pupil diameter was related to across-area latent variables.\nc-e, We next examined the changes in pupil diameter that occur just after saccade target onset (evoked or event-related pupil diameter). Here, we plotted the average raw pupil diameter trace from an example session from each animal (Monkey Pe: panel c, Monkey Wa: panel d, Monkey Sa: panel e). After target onset (time = 0), the diameter of the pupil fluctuated rapidly over the course of a few hundred milliseconds. To compensate for slow changes in pupil diameter from moment to moment, we identified a baseline period from 100-200ms before target onset (gray shading). The average baseline pupil diameter during this window was subtracted from the pupil diameter trace on a trial-by-trial basis. This subtraction equalized the size of the pupil at the start of each trial to effectively measure the event-related pupil diameter magnitude as a deviation from this baseline value. Note that this differs from the trial-by-trial pupil diameter measure, for which the baseline is identified on a moving average of 25 trials rather than on every trial.\nf-h, The baseline-subtracted pupil diameter traces averaged across trials are shown for an example session from each animal (Monkey Pe: panel f, Monkey Wa: panel g, Monkey Sa: panel h). We defined the event-related pupil diameter as the change in pupil diameter from the baseline period (-200 to -100ms before target onset, now at zero (a.u.)) to the manually identified event-related peak in pupil diameter (marked by arrows). This peak was manually adjusted based on the trial-averaged pupil trace in each session to account for shifts in observed pupil latencies across sessions (Monkey Pe: 221-297ms, Monkey Wa: 276-356ms, Monkey Sa: 127-200ms after target onset). We measured the mean pupil diameter in a 100ms window centered on the identified peak, which is relative to the baseline window on each trial. To compare across sessions, we z-scored the event-related pupil diameter within each session.\ni, We found that the magnitude of the event-related pupil diameter was negatively correlated with the trial-by-trial pupil diameter across sessions (mean Pearson's correlation, Monkey Pe: rho=-0.21, Monkey Wa: rho=-0.021, Monkey Sa: rho=-0.078) in two of three animals (p<0.0001; pooled across animals, n=41 sessions; Monkey Pe, p<0.0001, n=16; Monkey Wa, p=0.26, n=9; Monkey Sa, p=0.00026, n=16; one-tailed one-sample t-tests against zero). While the negative relationship we observed might appear counterintuitive, it in fact aligns with previous studies [8,14-16] and may reflect physical limits imposed by the musculature of the iris [17]. That is, during a time period that the trial-by-trial pupil diameter is especially large the biomechanics of the iris may constrain any event-related pupillary response to be smaller than usual (and vice versa).\nj, Although the trial-by-trial pupil diameter and event-related pupil might be correlated merely due to physical constraints from the iris musculature, it is also possible that they are influenced by different arousal processes and timecourses. As such, they might relate independently to the across-area component of neural activity. To this end, we performed an analogous regression analysis to that of Fig. 6 to relate the latent variables identified by pCCA-FA to the event-related pupil size (see Methods). Briefly, we regressed the full set of either the across- or within-area latent variables onto the event-related pupil diameter for each trial. We then assessed r2, the proportion of variance in event-related pupil diameter explained by the regression. As in Fig. 6c (same plotting conventions used here), we found that across-area latent variables explained more variance in the event-related pupil diameter than the within-area latent variables. We also report a 95% confidence interval over a “shuffle” distribution (gray bars, see Methods).\nAggregated across sessions, r2 values for explaining event-related pupil diameter aggregated across sessions were greater for across-area latent variables (mean r2= 0.035) than for within-area latent variables (right PFC, mean r2 = 0.013; left PFC, mean r2 = 0.013) for two of three animals (p<0.0001; pooled across animals and brain areas, n=82 brain areas, 2 per session; Monkey Pe: p<0.0001, n=32; Monkey Wa: p=0.43, n=18; Monkey Sa: p<0.0001, n=32; one-tailed t-tests on the differences in r2 for each session). When controlling for the number of latent variables (see Methods), we found similar results (p=0.00052; pooled across animals and brain areas, n=82 brain areas, 2 per session; one-tailed t-tests on the differences in r2 for each session).\nFurther, we determined that the across-area latent variables explained variance in the event-related pupil beyond what was explained by its correlation with the trial-by-trial pupil diameter (mean residual r2 = 0.028). To do this, we regressed out the contribution of the trial-by-trial pupil diameter to the event-related pupil, yielding a residual event-related pupil measure. Then, we assessed the variance explained by the across- and within-area latent variables in the residual event-related pupil diameter using the same linear regression approach as above (across-area residual r2 > shuffle r2, p<0.0001; pooled across animals, n=41 sessions; Monkey Pe, p<0.0001, n=16; Monkey Wa, p=0.13, n=9; Monkey Sa, p=0.0025, n=16; one-tailed one-sample t-tests against the mean of the shuffle distributions). Together with Fig. 6, we found that the across-area latent variables were related to multiple timescales of pupil activity, suggesting a link to one or more brainwide arousal processes.”\nMinor:\nL1360: “in a way” instead of “in way”\nFixed, thank you for catching this.\nL253: I am curious about the use of the word “complexity” here. Is it being used interchangeably with dimensionality? If so, is there a reason for that?\nWe thank the reviewer for this pointing this out. We use the term complexity as a more accessible way to refer to the shared dimensionality (d_shared) of the population, which we define on L224 (also used in Williamson, et al., 2016 and Bittner, et al., 2017). To clarify this definition, we have added a sentence to convey the more accessible terminology. However, we do note one place in the Methods section “Pupil diameter regression” where the terminology was inconsistent, so we replaced the term “dimensionality” with “number of latent variables”.\nChanges:\n● L226 in Results for Fig. 3:\n○ Addition: “A larger d_shared would indicate complex interactions among neurons, whereas a smaller d_shared would indicate less complex interactions among neurons.”\n● L952 in Methods section “Pupil diameter regression”:\n○ Original (in Methods): “To control for differences in across- and within-area dimensionality, we repeated the entire procedure using only one predictor (i.e., the top across- or within-area latent variable) rather than the full set of predictors (i.e., d across- or d_m within-area latent variables). This yielded similar results.”\n○ Original (from Fig. 6c caption): “We repeated the same analysis but used only the top co-fluctuation pattern of each type for the regression model to control for differences in across- and within-area dimensionality”\n○ Updated (combined in Methods): “To control for differences in the number of across- and within-area latent variables, we repeated the entire procedure using only one predictor (i.e., the top across- or within-area latent variable) rather than the full set of predictors (i.e., d across- or d_m within-area latent variables). This yielded similar results (p<0.0001; pooled across animals and brain areas, n=82 brain areas, 2 per session; one-tailed t-test on the differences in r2 for each session).”\nReviewer #2 (Remarks to the Author):\nThis paper develops a method combining canonical correlation analysis and factor analysis and applies it to simultaneously recorded spike activity from both sides of the monkey's lateral prefrontal cortex during a delayed eye movement task, and decomposes the intertrial variability of the activity into a common component between the left and right prefrontal areas and a component within each area. By showing that the common inter-areal component corresponds well to the pupil diameter reflecting the arousal level, they demonstrated that this method can explain specific phenomena effectively.\nBoth canonical correlation analysis and factor analysis are major data analysis methods, but I appreciate the fact that this method integrates them in a straightforward manner to obtain previously unobtainable results. The description is also clear and easy to understand, and the validity of the method and results is well verified. Therefore, I have no major objections to the scientific significance of this paper, and I will limit myself to making a few not-so-serious points regarding the acceptance or rejection of this paper. However, I have added some comments at the end that are not related to the acceptance or rejection of the paper. I hope they will find them useful in their future research.\nWe thank the reviewer for their positive comments, and appreciate their thoughtful suggestions for future research.\nMajor\nThe limitations of this method have not been discussed. What is supposed to be assumed by this method? What kinds of things might be missed? It seems natural, in a sense, that left and right hemisphere activity would co-vary with arousal state, but why was this not seen in a simple correlation and could only be obtained with the present method? etc. These questions need to be answered.\nWe thank the reviewer for these suggestions, we have reworked the Discussion in order to incorporate these answers.\nFirst, for the assumptions and limitations of pCCA-FA:\n1. Like other similar dimensionality reduction methods, pCCA-FA defines a linear-Gaussian relationship between the latent variables and the neural activity. We state this in the Methods section (L694), and we have now additionally brought this up as a potential extension of pCCA-FA in Discussion (L459).\n2. pCCA-FA is designed to identify latent variables that capture variance that is shared among neurons. Each method has a different objective function of what the latent variables are intended to capture. For example, one can use a targeted dimensionality reduction approach to identify activity that covaries with stimulus or behavioral variables (dPCA, Kobak, et al. 2016). The extension of pCCA-FA to identify other types of variance is now brought up in Discussion (L456).\n3. Limitations such as hyperparameter selection, run time, and other logistic concerns are now described more thoroughly in Discussion.\nSecond, for what could be missed:\n1. Related to the above assumption that pCCA-FA is a linear method, we could miss out in identifying nonlinear interactions among neurons. We have now noted this as an extension of pCCA-FA in Discussion (L459).\n2. pCCA-FA is a static method and, thus, does not incorporate the concept of time. pCCA-FA is designed to be applied to a single time bin of neural activity per trial (across many trials). Both of these limitations are not specific to pCCA-FA, but rather apply to the entire class of static linear-Gaussian dimensionality reduction techniques (including FA, pCCA, and principal component analysis). We have now noted this as an extension of pCCA-FA in Discussion (L452).\n3. Variance assigned to the within-area component of neural activity could, in fact, be shared with a non-recorded brain area. Similar interpretational considerations are present for any dimensionality reduction method. We bring this point up in Discussion (L385).\nThird, for why we could not uncover the connection to arousal using a simple correlation: A simple correlation measure (r ) was also able to uncover co-fluctuations across the hemispheres of PFC (Fig. 2). However, simple correlation does not isolate the interactions that are shared across hemispheres from those within hemispheres. By applying pCCA-FA, we specifically isolate the component of neural activity that is shared across hemispheres, and we found that this component (and not the within-area component) was related to arousal state (Fig. 6). As such, a simple correlation of neural activity to the pupil diameter would not have been as effective as our approach. We further discuss the limitations of mean r in Results section “Mean spike count correlation does not adequately partition shared variance”.\nChanges:\n● Reworked Discussion paragraph (L447):\n“pCCA-FA fills an unmet need in the set of exploratory tools tailored for studying the trial-to-trial variability of neural activity. When studying neurons from a single brain area, FA is among the simplest approaches for separating variance that is shared among neurons from variance that is independent to each neuron. For two brain areas, pCCA-FA is among the simplest methods for separating variance that is shared across areas from variance that is only shared among neurons within each area. To incorporate the concept of time, one can consider using Delayed Latents Across Groups (DLAG) [123] or related methods [33,124-129] to separate activity time courses that are shared across areas from those that are shared only among neurons within each area. pCCA-FA can be extended in the following ways. First, one may seek to identify latent variables that are explicitly related to task stimuli or the animal’s behavior by using a targeted dimensionality reduction approach [130-133] or by applying pCCA-FA to raw spike counts. Second, one may seek to identify non-linear interactions between brain areas [134-137]. As multi-area population recordings continue to increase in prevalence [28-30], pCCA-FA and related methods will become increasingly essential for dissecting the rich structure of across- and within-area interactions.”\nCanonical correlation analysis assumes that correlations exist. Will this method ever produce a pseudo correlation? In other words, will data from different trials or from different animals, measured under the same task conditions, produce any pseudo-positive results?\nWe agree with the reviewer that canonical correlation analysis and correlational methods as a whole are susceptible to pseudo-correlations. This motivated our preprocessing step of removing slow-timescale trends, as we illustrate in Supplementary Fig. 1. In Supplementary Fig. 1c, we show that it is indeed possible to uncover pseudo-correlations using pCCA-FA. This is not a drawback of pCCA-FA per se, but rather it applies to any correlational method. Thus, we must be cautious in our interpretation of correlations and perform the necessary controls and/or data preprocessing.\nIs it possible to interpret the significance of the dimensions obtained? If this question is wrong, there is no need to answer it.\nYes, this is an important question. To ensure that the co-fluctuation patterns identified by pCCA-FA are meaningful, we used cross-validation to determine the number of latent variables (i.e., the number of co-fluctuation patterns) present in the recorded activity (described in Methods). Briefly, we divide the trials from each session into a training set and a testing set, and train a pCCA-FA model at each of a wide range of candidate dimensionalities (1, 2, …,15 dimensions for both across- and within-area dimensionality in this work). Using the testing trials, we identified the candidate model (i.e., number of dimensions) that maximized the data likelihood. As such, we ensure that the dimensions identified in the pCCA-FA model generalized to the unseen testing trials, and are thus meaningfully capturing shared variance in the neural recordings. Further, our analysis linking the latent variables identified by pCCA-FA to the pupil diameter demonstrates that meaningful interactions in the brain were captured by the model. The latent variables are identified without any information about the pupil diameter. Thus, any relationship between the latent variables and the pupil diameter emerged naturally, rather than being “baked in” to the model.\nMinor\nI suppose they have performed single-unit sorting. If my understanding is correct, describe it clearly.\nWe have not performed single unit sorting. Rather, multi-unit activity is measured from each electrode on the Utah array. We report our procedure for removing waveforms not likely to be of neuronal origin in Methods. Thanks to the reviewer’s comment, we realized that single versus multi-unit differentiation was not mentioned, so we have updated the Methods section accordingly.\nChanges:\n● L536 in Methods section “Preprocessing of neural activity”:\n○ Addition: “The spikes from each electrode on the Utah array were treated as multi-unit activity, and no further differentiation of single units was performed. Thus, each electrode was treated as one neural unit.”\nIf rsc is calculated from multi-unit activity, will inter-regional correlations be obtained?\nPooling neurons on an electrode (i.e., multi-unit activity) tends to increase the magnitude of mean r (see Fig. 4 from Cohen & Kohn, 2011). This is because combining several units essentially averages out variability which is independent to each unit, yielding a large aggregate correlation between multi-units. In this work, we computed r from multi-unit activity for both across- and within-area pairs of electrodes. Pooling of neurons influences the across-and within-area r in a similar way, and thus it is fair to compare across- and within-area r that are both computed using multi-unit activity. We have added a statement to this effect in Methods.\nChanges:\n● L539 in Methods section “Preprocessing of neural activity”:\n○ Addition: “Using multi-unit activity tends to yield higher correlation values than single unit activity [8], and this applies equally to the across-area and within-area measurements.”\nIs it possible to interpret the significance of the dimensions obtained? If this question is wrong, there is no need to answer it.\nSee above response, this comment was repeated in the “Major” section.\nSchwartz and Goldman-Rakic's anatomical study [82] has indicated that nerve connections between the left and right prefrontal cortex are “comb-like.” The Utah array may also include areas that are directly innervated from the contralateral side and areas that are not. The correlations may then differ significantly depending on the electrodes in the Utah array.\nWe agree that the pattern of connections is important, and Schwartz & Goldman-Rakic’s study (now reference [96]) informs the PFC connectivity across hemispheres. We examined whether the co-fluctuation patterns identified by pCCA-FA exhibited any spatial clustering. For example, if patterns of across-hemisphere projections appeared on the array, we might expect to see clusters of electrodes with similar weights or periodic structure in the values of the weights. However, we did not find an apparent spatial structure within the scale of our array, which covers 4mm^2 of cortex in each hemisphere. Below, we provide a heatmap of the weights of the top across- (pink) and within-area (blue) co-fluctuation patterns from an example session (Monkey Sa). Positive and negative weights were spatially distributed throughout the array (electrodes removed during preprocessing are noted here with 0 weight in white). To address the reviewer’s comment, we have adjusted the wording of our reference to this citation in the Discussion section.\nChanges:\n● L378 in Discussion:\n○ Addition: “Second, variance shared across hemispheres may arise from direct anatomical projections [95], which has been shown to be spatially patterned [96].”\nIt would be helpful to the reader if the recovery of the %sv and d_shared in Fig. A3 had a figure that clearly explains why this method requires less data than pCCA. If it is not easy to illustrate, it is not required.\nWe thank the reviewer for this suggestion, and have added a new panel to Supplementary Fig. 4. This panel illustrates how pCCA-FA has fewer parameters than pCCA. We hope this helps to visualize why pCCA-FA performs better than pCCA in the limited data regime.\nChanges:\n● Supplementary Fig. 4b:\n● Caption:\n“b, We illustrate how pCCA-FA (top) has fewer parameters than pCCA (bottom). Suppose that we record from 10 neurons in each of areas 1 and 2. Together, the 20 neurons have a joint spike count covariance matrix (left). The color of each square in the heatmap indicates the corresponding value in the covariance matrix. The color scale is indicated below the joint covariance matrix, and applies to all heatmaps. Precise numerical values are arbitrary and for illustration purposes only.\npCCA-FA approximates this joint spike count covariance matrix as a sum of a low-rank across-area component (top left, W_1 and W_2), a block diagonal within-area component where each block is low-rank (top middle, L_1 and L_2), and a diagonal independent neuron component (top right, Ψ_1 and Ψ_2). This results in a parameter count of (n_1 + n_2) x d) for the across-area component, (n_1 x d_1) for the within-area 1 component, (n_2 x d_2) for the within-area 2 component, and (n_1 + n_2) for the independent neuron component.\nIn contrast, pCCA approximates the joint spike count covariance matrix as the sum of a low-rank across-area component (bottom left, W_1 and W_2) and a block diagonal matrix where each block is full rank (bottom right, Ψ_1 and Ψ_2). The block diagonal matrix includes all variance that is not captured by the across-area component, including any variance shared among neurons solely within each area and independent neuron variances. This results in a parameter count of (n_1 + n_2) x d) for the across-area component and (n_1 x n_1) + (n_2 x n_2) for the block diagonal matrix.\nThe number of parameters for the across-area component is the same for pCCA-FA and pCCA. However, pCCA-FA can have a substantially smaller number of parameters than pCCA in the remaining components, especially if the number of within-area latent variables in pCCA-FA is small relative to the number of neurons in each population. This explains the practical benefits of pCCA-FA over pCCA in the limited data regime (see panel a, Fig. 3e-f).”\nComments\nThis method is interesting as a data analysis method and has a wide range of applications. For example, it may be used to determine the extent to which fluctuations in U.S. and Japanese stock prices are due to bilateral and domestic factors, respectively. Of course, the method could be applied to LFP (local field potential) as well as to spike activity. I have no doubt about the usefulness of this method. However, I did not find in this paper much willingness or a clear perspective to delve into the nature of how the brain interacts between the cerebral hemispheres.\nWe appreciate the reviewer’s positive comment on the usefulness of pCCA-FA. We agree that our Discussion section was lacking in this regard. To address this, we have substantially reworked the first few paragraphs of the Discussion section to provide more perspective on the circuit motifs that give rise to across-area shared variance as well as cellular and circuit mechanisms that may underlie these motifs.\nChanges:\n● Reworked Discussion paragraph (L365):\n“Across-area variance may arise from a variety of circuit motifs. First, it may come from circuits that provide common inputs to both hemispheres of PFC. One such pathway involves the encoding of brainwide cognitive states, such as arousal. Arousal state modulation propagates through brainwide signaling that originates in the locus coeruleus [90,91]. The locus coeruleus projects to both hemispheres of PFC [92,93] and influences pupil diameter [94], a common readout of arousal state. Correlations between neural activity and pupil diameter have been observed across the brain [36,43,54], supporting the notion that the across-area variance we observed in PFC may be in part due to common inputs originating in the locus coeruleus. Another possible source of common input is sensory or motor signals. Each hemisphere of PFC represents both visual fields [17], and common motor plans have been shown to be encoded across the brain [36,41,43]. Second, variance shared across hemispheres may arise from direct anatomical projections [95], which have been shown to be spatially patterned [96]. Interactions across the hemispheres are central to various brain functions such as working memory [48,49], motor control [41,42], sensory perception [40], and spatial attention [56,57].”\n● Reworked discussion paragraph (L399):\n“Underlying the across- and within-area circuit motifs are various cellular and circuit mechanisms, which have been implicated in the shaping of correlations among neurons both across and within brain areas. First, the structure of inhibition is known to shape correlations among neurons, where the influence of inhibition depends on its spatial scale [26,104] and timing relative to excitation [105-108]. Second, changes at the level of the synapse, including synaptic depression [109] and connection strength [110], can also regulate correlated variability. Third, neuronal correlations are influenced by neuromodulators such as acetylcholine [111-113] and norepinephrine [91], which together have wide-ranging modulatory effects including altering the strength of thalamocortical inputs [114] and regulating global brain state [115]. Finally, the strength of excitatory synaptic connectivity has been linked to the strength of correlation [116]. Anatomically, long-range horizontal connections within a brain area [117-119], across brain areas [96], and across the corpus callosum [95], are all weaker than local recurrent connectivity. However, the strength of the detected correlations need not be directly linked to the strength of the anatomical connection between hemispheres or among distant neurons [120]. For example, weak neuromodulatory inputs could drive strong correlations among neurons by influencing the much stronger local circuit mechanisms that exist in local circuitry. Our analyses using pCCA-FA, like any other assessment based solely on spiking activity, do not explicitly distinguish between these cellular and circuit mechanisms. Instead, we focused on the spiking variability yielded by these mechanisms and dissected the components of spiking variability that are shared across versus within brain areas.”\nIf they wanted to examine interhemispheric interactions using the ODR task, they should have included an anti-saccade task or some other experimental paradigm that required specific information to be exchanged between the two hemispheres.\nWe agree with the reviewer in that the ODR task is a relatively simple experimental paradigm. However, we would argue that this task still does require the exchange and coordination of information across hemispheres (perhaps especially for vertical saccades). While the animals fixate centrally, information from both sides of the visual field about the location of the saccade target is communicated to each hemisphere of prefrontal cortex, with a contralateral bias. In order to process, store, and plan an eye movement to the remembered location, neural activity in both hemispheres must be coordinated. However, the point is well taken that an anti-saccade task and other related paradigms could be an especially powerful way to further probe interhemispheric communication.\nThe bin width of 1 second does not seem appropriate either. The exchange of information between hemispheres necessary to accomplish a task should occur in a much shorter time span. I know that many neurons in the prefrontal cortex have low firing rates. If so, then the instantaneous firing rate or other similar index should be used.\nWe agree with the reviewer that it’s important to match the timescale of analysis to the relevant timescale of the task. Our approach used 1s of the delay period, which was meant to aggregate the interactions within and across brain areas observed during the memory delay period. Trial-to-trial variability manifests across many timescales, and this somewhat slower timescale is in fact the largest contributor to the statistical interactions among neurons, while the millisecond-to-millisecond interactions tend to be statistically less noticeable, even among nearby neurons (Smith & Kohn, 2008). This has led spike count bins of several hundred milliseconds to a second being used in many studies of trial-to-trial variability (e.g., Afshar, et al., 2011; Ecker, et al., 2014; Athalye, et al., 2017; Khanna, et al., 2019; Stringer, et al., 2019). We have now added text to the Methods to indicate this rationale.\nChanges:\n● L564 in Methods section “Preprocessing of neural activity”:\n○ Addition: “This one-second window aggregated the interactions across and within brain areas observed during the delay period on the trial-to-trial timescale.”\nWhat information is carried from the contralateral area and how is it carried? Is it carried in different rhythms depending on the situation or is there a specific time window? Does a variety of information arrive unselected or is it carefully chosen? Some diseases, such as patients with apraxia, have problems with the consistency of information processing between hemispheres. It is the mission of brain scientists to try to get to the essence of how functional differentiation and integration between hemispheres is achieved.\nWe agree with the reviewer that these are important questions. To better place our study in the context of the broader work on inter-hemispheric coordination, we have reworked the Discussion paragraphs pasted above in response to the first comment in this section. We also note later in the Discussion (paragraph starting on L423) that it appears PFC applies some selection to the neural activity that is shared across hemispheres. Previous work has termed this gating a “communication subspace”, where only select signals are passed between hemispheres (Semedo, et al., 2019). And, indeed, understanding how the brain’s hemispheres interact will hopefully lead to treatments of diseases that involve the consistency of information processing between hemispheres.\nReviewer #3 (Remarks to the Author):\nIn their manuscript “Interactions across hemispheres in prefrontal cortex reflect global cognitive processing”, McDonnell et al. propose a method termed pCCA-FA (probabilistic canonical correlation analysis in combination with factor analysis) to study interactions within and between the left and right dorsolateral prefrontal cortex (also termed regions or areas here). For each neuron, pCCA-FA estimates the proportion of trial-by-trial spike count variance attributable to self, within-region, and across-region sources. The authors also apply SVD to extract the main latent variables underlying these sources. They show that this approach offers advantages over previous methods (rSC, pCCA, or FA). They applied this method to study the populations of each prefrontal cortex in rhesus monkeys trained on a memory-guided saccade task. They find that the primary latent variable representing the across-area co-fluctuation pattern is a more effective predictor of pupil diameter than the corresponding within-area latent variable. Since pupil diameter is a marker of arousal, a brain-wide process, this finding supports that pCCA-FA can find truthful and insightful patterns and sources of variability in the brain.\nI thought the approach the authors propose to study neuronal interactions is elegant and simple. pCCA-FA combines two well-established methods and, as far as I understand, is solidly constructed. They also provide an adequate justification for the need for a novel analysis tool. However, aspects related to the application of the method to the dataset and the interpretation of the results could be further refined.\nWe thank the reviewer for their positive remarks. We have incorporated the reviewer’s helpful feedback related to application of the method and interpretation of the results, as detailed below.\nGeneral points:\n1. There is a discrepancy between the title and narrative of the manuscript. From the title, one would think that the main focus of the paper is the study of “global cognitive processing,” which is indeed addressed in Results section 2.5, where pupil diameter is analysed. However, the discussion recurs only briefly to this topic. I would like to see more efforts to relate the findings to cognition or, alternatively, a change in the title.\nWe thank the reviewer for this suggestion, and agree that the Discussion was previously lacking in this regard. In response, we have substantially reworked the first few paragraphs of the Discussion to better place our results in the context of what is known about the biological sources of shared variance across areas and within an area, as well as the role of shared variance in cognitive processes. Below we summarize the changes we have made which support the use of the original title.\nChanges:\n● We substantially reworked the first part of the Discussion. We now devote a paragraph on the circuit motifs that may be a part of the across-area variance we observed. This paragraph also more extensively reviews the role of across-area variance in additional cognitive processes in the literature. We also devote a paragraph to cellular and circuit mechanisms that may underlie the across-area circuit motifs.\nL365: “Across-area variance may arise from a variety of circuit motifs. First, it may come from circuits that provide common inputs to both hemispheres of PFC. One such pathway involves the encoding of brainwide cognitive states, such as arousal. Arousal state modulation propagates through brainwide signaling that originates in the locus coeruleus [90,91]. The locus coeruleus projects to both hemispheres of PFC [92,93] and influences pupil diameter [94], a common readout of arousal state. Correlations between neural activity and pupil diameter have been observed across the brain [36,43,54], supporting the notion that the across-area variance we observed in PFC may be in part due to common inputs originating in the locus coeruleus. Another possible source of common input is sensory or motor signals. Each hemisphere of PFC represents both visual fields [17], and common motor plans have been shown to be encoded across the brain [36,41,43]. Second, variance shared across hemispheres may arise from direct anatomical projections [95], which have been shown to be spatially patterned [96]. Interactions across the hemispheres are central to various brain functions such as working memory [48,49], motor control [41,42], sensory perception [40], and spatial attention [56,57].”\nL399: “Underlying the across- and within-area circuit motifs are various cellular and circuit mechanisms, which have been implicated in the shaping of correlations among neurons both across and within brain areas. First, the structure of inhibition is known to shape correlations among neurons, where the influence of inhibition depends on its spatial scale [26,104] and timing relative to excitation [105-108]. Second, changes at the level of the synapse, including synaptic depression [109] and connection strength [110], can also regulate correlated variability. Third, neuronal correlations are influenced by neuromodulators such as acetylcholine [111-113] and norepinephrine [91], which together have wide-ranging modulatory effects including altering the strength of thalamocortical inputs [114] and regulating global brain state [115]. Finally, the strength of excitatory synaptic connectivity has been linked to the strength of correlation [116]. Anatomically, long-range horizontal connections within a brain area [117-119], across brain areas [96], and across the corpus callosum [95], are all weaker than local recurrent connectivity. However, the strength of the detected correlations need not be directly linked to the strength of the anatomical connection between hemispheres or among distant neurons [120]. For example, weak neuromodulatory inputs could drive strong correlations among neurons by influencing the much stronger local circuit mechanisms that exist in local circuitry. Our analyses using pCCA-FA, like any other assessment based solely on spiking activity, do not explicitly distinguish between these cellular and circuit mechanisms. Instead, we focused on the spiking variability yielded by these mechanisms and dissected the components of spiking variability that are shared across versus within brain areas.”\n● We reworked the motivation of the pupil analyses in Fig. 6 to provide background on the cognitive modulation of the pupil.\nL318: “Given that we found substantial across-area variance (Fig. 4) involving different co-fluctuation patterns than the within-area component (Fig. 5), we considered whether the across-area variance might be related to a brainwide process. Arousal is known to modulate activity of neurons throughout the brain [53,54,76-78]. In addition, arousal state has been linked to changes in both the tonic diameter of the pupil [79,80] and event-related changes in pupil diameter [81,82]. Changes in pupil diameter have been linked to arousal-related changes in neural activity in both cortical [83-86] and subcortical [79,87-89] regions. Given that arousal has been found to be a brainwide phenomenon, we wondered if it was related to the across-area variance we observed. Recall that the across-area component is defined as variance that is shared with one or more neurons in the other area (regardless of whether that variance is also shared with neurons in the same area). We therefore asked how well the across- and within-area latent variables identified by pCCA-FA explained variance in the pupil diameter (Fig. 6a).”\n● We added a new supplementary figure (Supplementary Fig. 8) which provides a more thorough description of event-related (evoked) pupil diameter and the cognitive modulation of this faster timescale of pupil activity.\n● Caption:\n“In Fig. 6, we investigated whether the across- or within-area components of neural activity were related to the trial-by-trial pupil diameter measured during the delay period (relative to the ongoing 25-trial average, see Methods). This assessed trial-to-trial fluctuations in pupil diameter relative to the ongoing mean (with a baseline subtraction that was matched to our neural data processing), and did not explicitly assess event-related changes in the pupil locked to saccade target onset. Measures of tonic pupil size have been linked to arousal and other cognitive states [10,11]. However, cognitive modulation has also been observed in evoked or event-related pupillary dynamics [12,13]. Here, we investigated whether the across- and within-area components of neural activity were related to event-related pupil diameter and compared that to the delay period pupil measurements presented in Fig. 6.\na, While animals performed the task, we monitored the raw pupil size during each trial's delay period over the course of the session (light gray). The size of the pupil slowly fluctuated over the course of the session [1], in this example session drifting upwards. As with our preprocessing of neural activity, we estimated the slowly varying component of pupil diameter by computing a moving average of 25 trials (dark gray, see Methods).\nb, We removed the estimated slow component of pupil diameter (top) from the raw pupil diameter. This yielded the residual “trial-by-trial” pupil diameter, one value per trial (bottom). In Fig. 6, we found that the trial-by-trial pupil diameter was related to across-area latent variables.\nc-e, We next examined the changes in pupil diameter that occur just after saccade target onset (evoked or event-related pupil diameter). Here, we plotted the average raw pupil diameter trace from an example session from each animal (Monkey Pe: panel c, Monkey Wa: panel d, Monkey Sa: panel e). After target onset (time = 0), the diameter of the pupil fluctuated rapidly over the course of a few hundred milliseconds. To compensate for slow changes in pupil diameter from moment to moment, we identified a baseline period from 100-200ms before target onset (gray shading). The average baseline pupil diameter during this window was subtracted from the pupil diameter trace on a trial-by-trial basis. This subtraction equalized the size of the pupil at the start of each trial to effectively measure the event-related pupil diameter magnitude as a deviation from this baseline value. Note that this differs from the trial-by-trial pupil diameter measure, for which the baseline is identified on a moving average of 25 trials rather than on every trial.\nf-h, The baseline-subtracted pupil diameter traces averaged across trials are shown for an example session from each animal (Monkey Pe: panel f, Monkey Wa: panel g, Monkey Sa: panel h). We defined the event-related pupil diameter as the change in pupil diameter from the baseline period (-200 to -100ms before target onset, now at zero (a.u.)) to the manually identified event-related peak in pupil diameter (marked by arrows). This peak was manually adjusted based on the trial-averaged pupil trace in each session to account for shifts in observed pupil latencies across sessions (Monkey Pe: 221-297ms, Monkey Wa: 276-356ms, Monkey Sa: 127-200ms after target onset). We measured the mean pupil diameter in a 100ms window centered on the identified peak, which is relative to the baseline window on each trial. To compare across sessions, we z-scored the event-related pupil diameter within each session.\ni, We found that the magnitude of the event-related pupil diameter was negatively correlated with the trial-by-trial pupil diameter across sessions (mean Pearson's correlation, Monkey Pe: rho=-0.21, Monkey Wa: rho=-0.021, Monkey Sa: rho=-0.078) in two of three animals (p<0.0001; pooled across animals, n=41 sessions; Monkey Pe, p<0.0001, n=16; Monkey Wa, p=0.26, n=9; Monkey Sa, p=0.00026, n=16; one-tailed one-sample t-tests against zero). While the negative relationship we observed might appear counterintuitive, it in fact aligns with previous studies [8,14-16] and may reflect physical limits imposed by the musculature of the iris [17]. That is, during a time period that the trial-by-trial pupil diameter is especially large the biomechanics of the iris may constrain any event-related pupillary response to be smaller than usual (and vice versa).\nj, Although the trial-by-trial pupil diameter and event-related pupil might be correlated merely due to physical constraints from the iris musculature, it is also possible that they are influenced by different arousal processes and timecourses. As such, they might relate independently to the across-area component of neural activity. To this end, we performed an analogous regression analysis to that of Fig. 6 to relate the latent variables identified by pCCA-FA to the event-related pupil size (see Methods). Briefly, we regressed the full set of either the across- or within-area latent variables onto the event-related pupil diameter for each trial. We then assessed r2, the proportion of variance in event-related pupil diameter explained by the regression. As in Fig. 6c (same plotting conventions used here), we found that across-area latent variables explained more variance in the event-related pupil diameter than the within-area latent variables. We also report a 95% confidence interval over a “shuffle” distribution (gray bars, see Methods).\nAggregated across sessions, r2 values for explaining event-related pupil diameter aggregated across sessions were greater for across-area latent variables (mean r2= 0.035) than for within-area latent variables (right PFC, mean r2 = 0.013; left PFC, mean r2 = 0.013) for two of three animals (p<0.0001; pooled across animals and brain areas, n=82 brain areas, 2 per session; Monkey Pe: p<0.0001, n=32; Monkey Wa: p=0.43, n=18; Monkey Sa: p<0.0001, n=32; one-tailed t-tests on the differences in r2 for each session). When controlling for the number of latent variables (see Methods), we found similar results (p=0.00052; pooled across animals and brain areas, n=82 brain areas, 2 per session; one-tailed t-tests on the differences in r2 for each session).\nFurther, we determined that the across-area latent variables explained variance in the event-related pupil beyond what was explained by its correlation with the trial-by-trial pupil diameter (mean residual r2 = 0.028). To do this, we regressed out the contribution of the trial-by-trial pupil diameter to the event-related pupil, yielding a residual event-related pupil measure. Then, we assessed the variance explained by the across- and within-area latent variables in the residual event-related pupil diameter using the same linear regression approach as above (across-area residual r2 > shuffle r2, p<0.0001; pooled across animals, n=41 sessions; Monkey Pe, p<0.0001, n=16; Monkey Wa, p=0.13, n=9; Monkey Sa, p=0.0025, n=16; one-tailed one-sample t-tests against the mean of the shuffle distributions). Together with Fig. 6, we found that the across-area latent variables were related to multiple timescales of pupil activity, suggesting a link to one or more brainwide arousal processes.”\n2. It is unclear to me why both the “mean pupil diameter” and the “pupillary evoked response magnitude” are analysed. The differences between these two measurements, as well as their cognitive interpretation, are not discussed.\nWe thank the reviewer for this point, and agree that the differences between the event-related pupil diameter (or pupillary evoked response, previously Fig. 6d) and the trial-by-trial pupil diameter (Fig. 6c) were not fully explained. In this work, the event-related pupil diameter (or evoked response) and trial-by-trial pupil diameter are measured during different windows of each trial. The trial-by-trial pupil diameter is measured during the delay period of each trial, while the event-related pupil diameter is measured just after target presentation (manually selected on each session, as shown in Supplementary Fig. 8f-h). In addition, they have a different baseline correction, with the trial-by-trial pupil diameter centered on the average of a 25 trial bin (in the same way as the neural activity) and the event-related pupil corrected by subtracting a pre-stimulus fixation period on the same trial. In order to better motivate the concepts of event-related pupil diameter, we moved Fig. 6d to a new supplementary figure devoted to the pupil (Supplementary Fig. 8, see the new figure pasted in the comment above). There, we present additional detail on how the pupil diameter is preprocessed, how event-related pupil diameter is defined, and its relationship to the trial-by-trial pupil diameter (Fig. 6c). In reworking the event-related pupil diameter analyses, we also adjusted the preprocessing of the pupil diameter. We adjusted the window in which we measured the trial-by-trial pupil diameter to further differentiate it from the event-related pupil diameter. The trial-by-trial pupil diameter is now exclusively measured during the delay period of each trial, from target offset to saccade initiation cue. We also added a filter for pupil data quality, which excluded one session from our analyses of the pupil (Fig. 6 and Supplementary Fig. 8). We have updated the text and Methods to reflect these changes and better motivate the measurement of event-related pupil diameter.\nChanges:\n● L321 in Results for Fig. 6:\n○ Addition: “In addition, arousal state has been linked to changes in both the tonic diameter of the pupil [79,80] and event-related changes in pupil diameter [81,82].”\n● L343 in Results for Fig. 6:\n○ Original: “Next, we asked whether the across- and within-area latent variables could explain variance in the pupillary evoked response magnitude, defined as the fast-timescale change in pupil diameter just after visual stimulus presentation (see Methods). We found that the across-area latent variables also explained more variance in the pupillary evoked response than the within-area latent variables (right PFC: p= 0.0001, left PFC: p<0.0001; Fig. 6d).”\n○ Updated: “Furthermore, we found that the across-area latent variables explained variance in the event-related changes in pupil diameter observed just after target onset (Supplementary Fig. 8; p<0.0001, pooled across animals and brain areas; one-tailed t-test).”\n● L585 in Methods section “Preprocessing of pupil diameter”:\n○ Original: “We measured pupil diameter in the time window starting 200ms before target onset until the saccade initiation cue on every trial.”\n○ Updated: “On each trial, we measured the pupil diameter during the delay period, from target offset until the saccade initiation cue.”\n● L621 in Methods section “Preprocessing of pupil diameter”:\n○ Addition: “To ensure pupil data quality, we computed the variance across trials of the trial-by-trial pupil diameter. We then removed outlier sessions that exceeded a variance threshold of two times the interquartile range of the variance values across sessions. This was performed for each animal individually. This procedure excluded one session from Monkey Wa from analyses of the pupil (Fig. 6, Supplementary Fig. 8), resulting in n=9 sessions for analysis.”\n3. In the discussion, it is mentioned that anatomical constraints on network circuitry may affect the correlation between neurons [line 363]. However, the cited papers are about visual, not prefrontal cortex. I would suggest to include more appropriate references, e.g. Kritzer & Goldman-Rakic, 1995, Journal of Comparative Neurology; Pucak, Levitt, Lund, & Lewis, 1996, Journal of Comparative Neurology\nWe thank the reviewer for this suggestion and have updated the citations at this location accordingly.\n4. The results for the two different eye-pupil measures appear to differ slightly, particularly for monkey “Wa.” Moreover, monkey “Wa” shows overall different values for rsc, %sv, and d_shared. Is there any interpretation of these differences?\nThe r2 for the event-related pupil diameter regression is indeed substantially weaker than the r2 for the trial-by-trial pupil diameter regression for all animals, and especially for Monkey Wa. This could reflect different timescales of arousal-linked modulation of the pupil. It may be that brainwide arousal signals are most strongly linked to the pupil on the timescale of several seconds, the difference between the delay period pupil and the 25-trial average (the trial-by-trial pupil). The event-related pupil captures pupillary changes on a timescale of a few hundred milliseconds (the evoked response to the target relative to the fixation baseline), which reflects changes due to the luminance transient as well as arousal-linked modulation.\nMonkey Wa also tended to be an outlier in our other results such as r (Fig. 2d), %sv (Fig. 4b), and d_shared (Fig. 4c). We attribute this to two factors: 1) Monkey Wa tended to have the least number of trials in each experimental session (L527) which led to noisy estimation of latent variables and 2) a small number of neurons in one hemisphere (L555) which could limit the usefulness of estimating across-area interactions. We report all of the per-session results for all 3 animals and report p-values individually in each figure caption to fully convey this variation among animals.\nParticular points:\n[103-105] The phrasing in this sentence could be clearer. Since coordination between neurons is defined as synchronized patterns of activity, describing it as a process that causes or “results” in such patterns may appear tautological. This could also unintentionally suggest that neurons coordinate through signals other than their activity. Rewording the sentence could help avoid this ambiguity.\nWe thank the reviewer for pointing this out. We have updated the text at this location accordingly.\nChanges:\n● L104 in Results:\n○ Original: “The coordination among neurons during perception, cognition, and action results in aspects of each neuron’s activity being shared with neurons in the same brain area as well as in other brain areas.”\n○ Updated: “To support perception, cognition, and action, each neuron’s activity can be shared with other neurons in the same brain area, as well as with neurons in other brain areas.”\n[140] What statistical test was performed here? If I understood correctly, none of the tests shown in Fig. 2 directly compare the distributions of rSC within versus across regions.\nFrom the reviewer’s question we realized in hindsight that we could have been clearer about the statistical tests that were performed in each analysis. We actually did have a test that directly compared the across-area mean r to the within-area mean r in all sessions using paired one-tailed t-tests. However, the mention of the statistical test was only in the Fig. 2d caption and not the main text. To address this, we have updated our notation in each location a p-value is mentioned in the main text to include more details about the statistical test.\nChanges:\n● Example L141 in Results: (p-value notation was updated in many locations)\n○ Original: (Fig. 2d, p<0.0001)\n○ Updated: (Fig. 2d; p<0.0001, pooled across animals; paired one-tailed t-test)\nAdditionally, for monkey “Wa,” the distributions of average rSC within and across regions across all sessions do not appear very different to the naked eye.\nIndeed, when comparing the across-area mean r to the within-area mean r for Monkey Wa, there were some sessions in which the within-area r was not greater than the across-area mean r . This could reflect the sampling of neurons from each brain area, which was unbalanced in number for this animal. We note the p-values for comparing across-area to within-area mean r for each monkey in the caption for Fig. 2d. For Monkey Wa, we found that across all sessions, within-area mean r was greater than across-area mean r (p=0.039; pooled across brain areas, n=20 brain areas, 2 per session; one-tailed t-tests on the differences in mean r for each session).\n[248-251] Is this discrepancy the only reason that rSC is smaller than %sv? Could it be resolved simply by taking the absolute value or the squared value (as is done with %sv)?\nYes, because r is signed, positive and negative r values will effectively cancel out in a mathematical average. Because %sv incorporates positive and negative weights, it does not cancel out negative correlations among neurons. There is a mathematical relationship between r and the parameters of the pCCA-FA model (%sv). This relationship is analogous to the relationship derived in Math Note A of Umakantha, Morina, Cowley, et al., 2021 (reference [72]). We have added a citation to this reference in the results section on the limitations of r .\nChanges:\n● L261 in Results for Fig. 4:\n○ Original: “This discrepancy can be understood by considering the sign of pairwise correlations among neurons (Extended Data Fig. A2).”\n○ Updated: “This discrepancy can be understood by considering the sign of pairwise correlations among neurons (Fig. 2e, see Math Note A of [72]).”\n[316-317] The “pupillary evoked response magnitude” is measured “just after visual stimulus presentation,” which occurs before the delay period. From the current wording (“across-area latent variables were able to explain multiple aspects of pupil-linked arousal”), one might infer a causal role from the latent variable to the pupillary evoked response. However, the neuronal activity is recorded “1 second of the delay period immediately preceding the earliest possible saccade initiation cue in the session.” Therefore, it cannot be concluded that the “latent variables identified by pCCA-FA explained changes in pupil diameter,” since the change in pupil diameter precedes the neuronal activity, implying that the causal direction would be reversed.\nWe thank the reviewer for pointing out the discrepancy in the language surrounding the analysis of the pupil, and have updated it accordingly. We agree that we cannot infer causality from this analysis. In Fig. 6 and Supplementary Fig. 8j (previously Fig. 6d), we apply linear regression to uncover a correlational relationship between the latent variables identified by pCCA-FA and our two measures of pupil size. We then assessed r2, the statistical measure of variance explained by the regression models. We have clarified the wording, and now state that the latent variables explain variance in the pupil diameter from trial-to-trial.\nChanges:\n● L330 in Results for Fig. 6:\n○ Original: “We therefore asked how well the latent variables identified by pCCA-FA explained changes in pupil diameter (Fig. 6a).”\n○ Updated: “We therefore asked how well the latent variables identified by pCCA-FA explained variance in the pupil diameter (Fig. 6a).”\n● L343 in Results for Fig. 6:\n○ Original: “Next, we asked whether the across- and within-area latent variables could explain variance in the pupillary evoked response magnitude, defined as the fast-timescale change in pupil diameter just after visual stimulus presentation (see Methods). We found that the across-area latent variables also explained more variance in the pupillary evoked response than the within-area latent variables (right PFC: p= 0.0001, left PFC: p<0.0001; Fig. 6d). Thus, the across-area latent variables were able to explain multiple aspects of pupil-linked arousal.”\n○ Updated: “Furthermore, we found that the across-area latent variables explained variance in the event-related changes in pupil diameter observed just after target onset (Supplementary Fig. 8; p<0.0001, pooled across animals and brain areas; one-tailed t-test). Thus, the across-area latent variables may be related to multiple aspects of pupil-linked arousal.”\n[331-332] I would argue that this is not, in general, “the shared variance among PFC neurons.” Here, task conditions (right or left stimulus cue, or right or left response) are pooled together. Therefore, it specifically measures the task-condition-independent variance. This distinction is significant because the PFC is primarily involved in processes that are task-dependent, where the task condition plays an important role. Thus, the analysis is capturing the shared variance of a specific sub-process within the PFC.\nWe agree with the reviewer that this statement was an oversimplification of the results, and have updated the text accordingly.\nChanges:\n● L355 in Discussion section:\n○ Original: “Applying a novel statistical approach called pCCA-FA, we found that of the substantial shared variance among PFC neurons, a large proportion involved neurons in both areas.”\n○ Updated: “Applying a novel statistical approach called pCCA-FA, we found that of the substantial shared trial-to-trial variability among PFC neurons, a large proportion involved neurons in both areas.”\n[481-485] Why was the entire electrode signal removed instead of only removing the repeated units? Could this approach be removing nearby neurons that might be more connected, thereby further reducing the %sv of inter-region units?\nWe appreciate this suggestion. In this work, we analyzed multi-unit activity, which is the aggregate activity recorded on one electrode of the Utah array. As we did not sort for single units, it was not possible to remove a single unit off of an electrode. We agree that this was not clearly stated in the Methods section, and have updated it accordingly. It is possible that neglecting single unit activity reduces our sensitivity to measure %sv. However, as with any measurement technique, the spatial scale of the measurement determines the sensitivity of the result. We have now noted this in the results section.\nChanges:\n● L536 in Methods section “Preprocessing of neural activity”:\n○ Addition: “The spikes from each electrode on the Utah array were treated as multi-unit activity, and no further differentiation of single units was performed. Thus, each electrode was treated as one neural unit.”\n● L265 in Results for Fig. 4:\n○ Addition: \"Note that the values of across- and within-area %sv (like mean r ) can depend on neuronal properties, including their physical distance and tuning preferences”\n[500-502] If I understand correctly, the intention of this procedure is to remove condition-specific rate properties. However, this method only removes cue-specific rate variations, omitting response-specific variations. In other words, when pooling both correct and incorrect trials, for the same cue condition, there can be two overlapping modes of firing rate depending on whether the response was right or left. The paper does not specify which trials were selected (all or only correct) or what the hit rate was.\nWe agree. In this work, we analyzed only correct trials. We appreciate the reviewer pointing out that this was not specified, and we have added a sentence in the Methods section to address this. It is true however that variability in the precise motor action (though the saccade target location is identical) is always present in any experiment.\nChanges:\n● L525 in Methods section “Behavioral task”:\n○ Original: “The number of trials per session for Monkey Pe was 729.8 +/- 173.4 (mean +/- standard deviation), for Monkey Wa was 463.5 +/- 170.1, and for Monkey Sa was 850.1 +/- 284.2.”\n○ Updated: “For each session, we considered only trials in which the animal successfully acquired the saccade target (i.e., correct trials). The number of correct trials per session for Monkey Pe was 729.8 +/- 173.4 (mean +/- standard deviation), for Monkey Wa was 463.5 +/- 170.1, and for Monkey Sa was 850.1 +/- 284.2.”\nFurthermore, it would strengthen the manuscript to demonstrate that condition-mean subtraction effectively removes condition-specific activity patterns. This could be done, for example, by training a model to classify trials by condition before and after rate manipulation, showing that condition-mean subtraction renders the conditions indistinguishable.\nWe thank the reviewer for this suggestion. Our main purpose in combining across conditions (saccade target direction) was to boost the number of trials for model fitting (tens of trials per condition). If we had hundreds of trials per condition, this step would be unnecessary. We completed the suggested analysis by performing linear discriminant analysis in a cross-validated manner to predict the condition label before and after removal of condition-specific means. We indeed found that cross-validated decoding accuracy of the condition label before removal of condition-specific means (51.9% averaged across n=42 sessions) was much higher than the decoding accuracy after removal of condition-specific means (12.0% averaged across n=42 sessions). However, we note that the decoding accuracy after removing condition-specific means was still greater than chance (25%, 12.5%, or 6.25% depending on the session). This is likely due to the fact that, for a given neuron and condition, the distribution of spike counts is Poisson-like. Thus, the variance of the residuals after subtracting the condition-specific mean still can reflect the mean spike count. The reviewer makes a good point in that this could contribute to the shared variance identified by pCCA-FA (among neurons, shared across or within areas). We updated the text in the Methods section to note this. Additional details related to this point are included in response to the comment below.\nChanges:\n● L586 in Methods section “Combining trials across conditions”:\n○ Original: “To increase the number of observations for model fitting, we combined trials across conditions (saccade target direction). To accomplish this, before fitting the pCCA-FA model or before computing r , we first subtracted the condition mean from the binned spike counts in each condition and then combined residual spike counts across conditions.”\n○ Updated: “To increase the number of observations for model fitting, we combined trials across conditions (saccade target direction). To focus our analyses on trial-to-trial variability, for each neuron, we first subtracted the condition mean from the binned spike counts in each condition and then combined residual spike counts across conditions. Because the spike count variance can differ across conditions, the residual spike counts can still contain a small amount of condition-specific signal.”\n[503] It is mentioned that z-scoring “yielded qualitatively similar rSC results,” but nothing is said about its effect on the pCCA-FA results. Should we expect that z-scoring would also lead to similar outcomes in pCCA-FA as in rSC?\nWe appreciate the reviewer bringing this up, especially in conjunction with the previous comment. This has motivated us to repeat the key analysis related to the identification of shared variance among neurons (Fig. 4). In the main figure, we remove only condition averages before combining across conditions. Below, we also divide the spike counts of each neuron by its condition-specific variance (i.e., z-scored) before combining across conditions and fitting pCCA-FA. We found qualitatively similar results. As in Fig. 4, we found that across-area %sv and d_shared were greater than the within-area %sv and d_shared, respectively (%sv: p=0.0056, d_shared: p=0.00013; pooled across animals, n=84 brain areas, 2 per session; paired one-tailed t-tests). This is shown in the figure below, which is analogous to Fig. 4b-c in the manuscript.\nChanges:\n● L591 in Methods section “Combining trials across conditions”:\n○ Original: “An alternative approach is to z-score spike counts in each condition before combining, which yielded qualitatively similar r results”\n○ Updated: “An alternative approach is to z-score the spike counts for each neuron and condition before combining across conditions [63]. We repeated the key analyses related to the identification of shared variance among neurons in Figs. 2 and 4 using condition-specific z-scoring and found qualitatively similar results.”\n[818-827] Why was the null distribution of the regressions created by shuffling neuronal activity from neurons in other sessions? Why not simply shuffle the pupil diameter values within a session while using the firing rates of the same neuron but from other trials?\nWe appreciate this suggestion, and agree that it would also be a valid control. We have found that using time-matched trials from another session is a stricter null distribution than shuffling trials within sessions. Below, we repeated the analysis in Fig. 6c as a demonstration that shuffling the pupil values within each session yielded qualitatively similar results. The black plotted points are identical to those in Fig. 6c, and the gray bars indicating the 95% confidence interval of the shuffle distribution here is the within-session shuffle. We found that the null distributions in our original shuffle distribution were slightly wider than those plotted below, indicating that our original approach yielded a more conservative null distribution. However, we agree that either approach could be appropriate.\nChanges:\n● L948 in Methods section “Pupil diameter regression”:\n○ Addition: “An alternative method to generate the chance distribution involves performing a within-session shuffle, in which the latent variables from session i are used to explain variance in the pupil diameter from session i but with permuted trial labels. This procedure yielded a similar, but slightly more lenient, chance distribution.”\nReviewer #4 (Remarks to the Author):\nThis is an interesting and important paper of potentially high impact describing a new statistical technique developed by the authors that can parse components of variation in population activity patterns attributable to local interactions between neurons and inputs from neurons in another brain area.\nThis is a key challenge in deciphering the computational organization of the cerebral cortex, and it is the advantage of this technique over alternative methods. The technique could significantly advance our understanding of the interplay between local processing in cortical areas and large scale communication in cortical networks.\nNearly all of the questions I initially had about the validity of the technique were answered by results provided in Extended data. Generally, I recommend promoting as many of the highlights of the Extended results into additional panels of the main figures as possible, so that the majority of the readership sees the strong evidence validating the new technique.\nIn addition to describing and validating a new analytical technique, the paper provides important new information about the coordination of neural activity in prefrontal cortex, including the findings that correlations between brain areas (left and right hemispheres) are stronger than correlations within areas, and that the strength of within and between area neural correlations track several behavioral variables.\nWe thank the reviewer for their positive comments, and appreciate the suggestion to promote significant panels of the supplements to the main figures. We reworked both Figure 2 and 3 to this point, see below for details.\nMajor\n1. The paper could be strengthened by enhancing the presentation of the data validating the pCCA-FA technique. Much of the crucial information is provided in the Extended Data. This includes: (1) comparison of detected correlations in the original data with trial shuffled data (Fig. A1), (2) relation of noise correlation to signal correlation (Fig. A2), (3) recovery of ground truth parameters from simulated data (Fig. A3), and (4) evidence that the pCCA-FA technique outperforms alternative methods (Fig. A3). The paper could be strengthened if some of the key panels from these figures were promoted to the main figures of the paper (for example adding them to Fig. 3). This would strengthen the case that the analysis developed by the authors measures real neural correlations, successfully parses variation within and between cortical areas, and outperforms existing methods\nWe thank the reviewer for these suggestions to promote more results to the main text - we agree. The validation analyses were originally left to the supplements, as we tried to balance their inclusion in the main text with maintaining the flow of the main text. In response to the reviewer’s suggestions, we have made the following changes: promoting the signal correlation result to Fig. 2, and promoting one of our model validation simulations to Fig. 3 which includes a comparison to the alternative method pCCA. Furthermore, we have added a statement in the text that directly references the trial-shuffled result in Supplementary Fig. 1. We feel that the message in Supplementary Fig. 1 is strengthened by having the space to fully describe the procedure and results. As such, we opted to keep this figure in the supplement (please see the comment below where we address the trial shuffling analysis in more detail).\nChanges:\n● L258 in Results for Fig. 4:\n○ Original: “Thus, substantial across-area interactions were present in the activity and identified by pCCA-FA, contrary to the intuition from measuring pairwise across-area r .”\nsc\n○ Updated: “Thus, substantial across-area interactions were present in the PFC activity and identified by pCCA-FA (see also Supplementary Fig. 1e for the across-area correlation values computed using pCCA-FA), contrary to the intuition from measuring pairwise across-area r .”\nsc\n● Figure 2:\n○ Original: signal tuning result was presented in Supplementary Fig. 2.\n○ Updated: signal tuning result is now featured in Fig. 2e.\n● L150 in Results for Fig. 2:\n○ Addition: “Indeed, when we computed mean r among subsets of across-area neuron pairs with similar tuning (Fig. 2e, top) and opposite tuning (Fig. 2e, bottom), we found that similarly tuned pairs of neurons exhibited positive noise correlation (mean r = 0.019; p<0.0001, pooled across animals; one-tailed t-test), whereas oppositely tuned pairs of neurons exhibited negative noise correlation (mean r = −0.013; p<0.0001, pooled across animals; one-tailed t-test)”\nsc\n● L569 in Methods section “Computing signal correlation”:\n“In Fig. 2e, we investigated if the distribution of across-area r , despite having a mean near zero, contained pairs of neurons with positive r and other pairs of neurons with negative r , in a way that was systematically related to the neurons’ task tuning. To do this, we identified each neuron’s tuning across all of the saccade target locations. We determined whether each pair of neurons had similar or opposite tuning by computing their signal correlation. This was defined as the Pearson correlation between the two neurons’ raw trial-averaged responses to each of the possible targets (4, 8, or 16 depending on the session) during the one-second window described above. To assess significance for each neuron pair, we generated a chance distribution of correlation values by using a permutation test (shuffling the target location labels). We then labeled a neuron pair as having significant signal correlation (positive for similarly tuned and negative for oppositely tuned pairs) if their signal correlation was more extreme than the 99th percentile of the chance distribution. Neuron pairs not meeting this criterion were excluded from this analysis.”\n● Figure 3:\n○ Original: the validation of pCCA-FA and pCCA as a function of trial count was presented in Supplementary Fig. 3.\n○ Updated: this is now featured in Fig. 3e-f.\n● L892 in Methods section “Simulated data generation”:\n“For model validation in Fig. 3e-f, we defined these parameters in the following way. We generated random loading matrices Wm and Lm for each area m using independent draws from a standard normal distribution. We scaled the %sv in each loading matrix by adjusting the singular values of the across- and within-area loading matrices (W_m and L_m) relative to each other to achieve the desired ratios (~15% across-area and ~10% within-area). The independent variance values (i.e., entries of Ψ_m) were then scaled to achieve a level of 75%. Each simulated dataset had 30 neurons per area, 5 across-area latent variables (d = 5), and 3 within-area latent variables per area (d_m = 3 for each area m). For each trial count configuration, we generated 30 sets of ground truth parameters and a single simulated dataset per set of parameters.”\n● L915 in Methods section “Fitting pCCA-FA and pCCA to simulated data”:\n“To fit pCCA-FA to simulated data for the analysis in Fig. 3e-f, we used the EM algorithm described above. Across- and within-area dimensionalities were estimated in a cross-validated manner. Candidate dimensionalities for across- and within-area ranged from 1 to 10 (inclusive), for both across- and within-area latent variables. After the across- and within-area dimensionalities were selected, a final pCCA-FA model was fit with the specified dimensionalities using all simulated trials. The parameters of this final model were used to compute %sv and d_shared.\nTo fit pCCA, we used an analogous EM algorithm, presented in [74]. We fit pCCA to the same sets of simulated data as pCCA-FA. The across-area dimensionality was estimated in a cross-validated manner. Candidate dimensionalities for the across-area latent variables ranged from 1 to 10 (inclusive). After the across-area dimensionality was selected, a final pCCA model was fit with the specified dimensionality using all simulated trials. The parameters of this final model were used to compute across-area %sv and d_shared.”\n2. Related to the above, trial shuffling is utilized in the extended data to validate removal of slow sources of correlation over trials (Fig. A1). The current data could be strengthened by showing that trial shuffling destroyed detected neural correlations at faster time scales after removal of the slower components of correlation. This would substantiate that detected neural correlations depended on the simultaneity of the neural signals and were not artifacts of the analysis. The amount of correlation detected in the trial shuffled data (in some cases ~10%, Fig. 1Ad) was not trivial in relation to the magnitude of the correlation in the original data (~10-15%, Fig. 4B). It would therefore be important to show that the neural correlations detected by this method significantly exceed the level of correlation seen after trial-shuffling.\nWe agree that the current data can be strengthened in this way. Following the reviewer’s suggestion, we computed the cross-validated canonical correlations from the fast component of non-shuffled neural activity. We found that the cross-validated canonical correlations in this case (non-flipped) were all higher than those shown in Supplementary Fig. 1d (flipped). As the reviewer points out, this indicates that the across-area correlations detected by pCCA-FA depend on the simultaneity of the neural recordings. We have added this new analysis to Supplementary Fig. 1 (panel e).\nChanges:\n● L258 in Results for Fig. 4:\n○ Original: “Thus, substantial across-area interactions were present in the activity and identified by pCCA-FA, contrary to the intuition from measuring pairwise across-area r .”\nsc\n○ Updated: “Thus, substantial across-area interactions were present in the PFC activity and identified by pCCA-FA (see also Supplementary Fig. 1 for the across-area correlation values computed using pCCA-FA), contrary to the intuition from measuring pairwise across-area r .”\nsc\n● Supplementary Fig. 1e:\n● Caption:\n“Throughout this work, we fit pCCA-FA to the residual fast component of neural activity. Here, we sought to assess whether the cross-validated canonical correlations (ρ fast) in this case were meaningful. To do this, we defined a chance distribution using the canonical correlations recovered when fitting pCCA-FA to the flipped fast component (ρ fast-flip, i.e., the values plotted in panel d). We found that the cross-validated canonical correlations in the fast component of our PFC recordings were significantly above chance (ρ fast − ρ fast-flip > 0) across all sessions (Monkey Pe: p < 0.0001, n = 103 co-fluctuation patterns; Monkey Wa: p < 0.0001, n = 28; Monkey Sa: p < 0.0001, n = 132; paired one-tailed t-tests). This indicates that substantial trial-to-trial correlation exists between the hemispheres of PFC.”\n3. A key claim/finding is that the new analysis can parse correlations driven by local interactions between neurons from correlations driven by input from other brain areas. If this is true, differences in %sv attributable to interactions within and between cortical areas should be greater in the original data than after shuffling neurons between brain areas (randomly changing right/left hemisphere labels). Demonstrating that with an additional bootstrap test if feasible would strengthen the main claim of the paper.\nWe thank the reviewer for this suggestion. We have run the suggested control analysis where we shuffled the hemisphere labels of our neural recordings and fit a pCCA-FA model to these shuffled spike counts. In this case, we would expect within-area shared variance to be significantly reduced. Within-area variance is only shared among neurons within one brain area, and not with any neurons in the other area. Shuffling the labels of the neurons in each area would destroy this shared variance. Indeed, the amount of within-area %sv was reduced to near zero in the shuffled data. The across-area %sv remained comparable to its original level. This is because across-area %sv is shared among neurons in both areas, so the original hemisphere label is irrelevant. We have added this analysis as a new Supplementary Fig. 5 (pasted below, right panel). For comparison, we have also included here a copy of Fig. 4b (pasted below, left panel), which is the unshuffled (i.e., true hemisphere labels) analog of the new Supplementary Fig. 5.\nChanges:\n● Supplementary Fig. 5 (above, right panel)\n● Caption:\n“In Fig. 4b, we compared the amount of shared variance (%sv) in the across-area component versus the within-area component. Here, we validated that we identified meaningful shared variance in our PFC recordings by randomly permuting the area labels (i.e., right or left hemisphere) of the neurons. In other words, we randomly assigned each neuron to area 1 or area 2. When we permute the area labels of the neurons, each shuffled population contains some true area 1 neurons and some true area 2 neurons. We can think of each shuffled population as comprising neurons drawn from the same overall population of neurons. Thus, it is unlikely for there to be activity variance that is shared exclusively among neurons within a shuffled population. We would expect nearly all of the shared variance among neurons to be assigned to the \"across-area\" component, and little shared variance to be assigned to the \"within-area\" component.\nTo perform the shuffle control, we took the same spike counts that were used to fit pCCA-FA (see Methods) and permuted the area labels of the neurons. Then, we fit pCCA-FA to the spike counts of the shuffled neurons in a cross-validated manner, identical to the analyses in the main figures (Figs. 4-6, see Methods). We computed the across- and within-area %sv using the parameters of the resulting pCCA-FA model in each session. Compared to Fig. 4b (which is the same analysis using the pCCA-FA models fit using true area labels), the within-area %sv (horizontal axis) in the shuffled models was reduced to near zero. The across-area %sv (vertical axis) was similar or slightly greater for the shuffled neurons compared to the neurons with true area labels. This indicates that the true area labels of the neurons and the differences we found in the across-area versus within-area interactions in Fig. 4 are meaningful.”\n4. More %sv was attributed to interactions between than within brain areas. The models derived to estimate %sv between areas had more free parameters (a greater number of latent parameters). It will be important for the main findings of the paper to argue/show that the difference in %sv attributed to local and long tract interactions did not reflect differences in the number of latent parameters in the models used to fit the two types of interactions.\nWe agree with the reviewer that showing this is essential to the claims of our work. In Supplementary Fig. 2, we demonstrate in simulation that regardless of whether the across-area or within-area interactions dominate in the ground truth model, pCCA-FA is able to accurately identify the proportion of variance attributed to the across- and within-area components. For example, consider the analysis in Supplementary Fig. 2a, where the number of across-area latent variables (d=5, indicating a larger number of free parameters) is higher than the number of within-area latent variables (d_m=3, indicating a smaller number of free parameters). If there were a bias in estimation due to the number of parameters, we might find that across-area %sv is over-estimated and within-area %sv is under-estimated. Instead, in the lower panel of Supplementary Fig. 2a, we show that we can identify %sv accurately when across-area %sv is greater than (above dashed line), equal to (on dashed line), or less than (below dashed line) the within-area %sv. This indicates that the number of latent variables, which was greater for the across-area component in all three cases shown in Supp Fig. 2a, does not bias our ability to accurately identify the level of %sv. To clarify this in the manuscript, we have added a statement about this control in the Results when describing Fig. 4.\nChanges:\n● L273 in Results for Fig. 4:\n○ Addition: “We further note that these findings are not a consequence of statistical bias in measuring across- versus within-area %sv and d_shared (Supplementary Fig. 2).”\n5. The analysis is built to detect noise correlations, e.g. correlations in neural signals after modulation of activity due to task performance has been subtracted out. Relations to behavior (pupil diameter, spatial tuning) are therefore indirect. Can the technique be adapted to detect correlations in neural signal (task-related modulations in firing rate)? This would enhance links between detected correlations and behavior. Some description of how this might be done in future studies in the Discussion could be useful to understand the potential applicability of the technique to different questions.\nYes, the method can be applied to raw spike counts (e.g., task-related modulations in firing rate) as well as trial-to-trial variability of spike counts for a given condition (which we used in this work). If we applied pCCA-FA to raw spike counts in PFC, we might expect to find saccade target information in both the across- and within-area components. The proportion of target-related activity in each component might depend on which hemifield the target was in. This approach could provide further insight into inter-hemispheric interactions in future studies. We have updated the Discussion to clarify the application of pCCA-FA to other types of neural activity.\nChanges:\n● Reworked Discussion (L455):\n“pCCA-FA can be extended in the following ways. First, one may seek to identify latent variables that are explicitly related to task stimuli or the animal’s behavior by using a targeted dimensionality reduction approach [130-133] or by applying pCCA-FA to raw spike counts. Second, one may seek to identify non-linear interactions between brain areas [134-137]. As multi-area population recordings continue to increase in prevalence [28-30], pCCA-FA and related methods will become increasingly essential for dissecting the rich structure of across- and within-area interactions”\n6. It would be useful to clarify the link between detected neural correlations and the physiological processes that are likely to generate them, namely synaptic interactions between neurons. In particular, how do the authors envision that synaptic interactions between neurons could drive correlations in firing rate at the trial level that are captured by their analysis? For example, correlations in delay period firing rate over trials could reflect net excitatory and inhibitory synaptic interactions between neurons integrated over the 1 sec delay period, etc.\nWe thank the reviewer for this suggestion. We agree that the synaptic interactions among neurons likely underlie the measured neuronal correlations, and were not adequately discussed. To address this, reworked the first few paragraphs of the Discussion, which now covers the various circuit motifs and cellular and circuit mechanisms that have been shown to shape correlations among neurons.\nChanges:\n● Reworked Discussion paragraph (L365):\n“Across-area variance may arise from a variety of circuit motifs. First, it may come from circuits that provide common inputs to both hemispheres of PFC. One such pathway involves the encoding of brainwide cognitive states, such as arousal. Arousal state modulation propagates through brainwide signaling that originates in the locus coeruleus [90,91]. The locus coeruleus projects to both hemispheres of PFC [92,93] and influences pupil diameter [94], a common readout of arousal state. Correlations between neural activity and pupil diameter have been observed across the brain [36,43,54], supporting the notion that the across-area variance we observed in PFC may be in part due to common inputs originating in the locus coeruleus. Another possible source of common input is sensory or motor signals. Each hemisphere of PFC represents both visual fields [17], and common motor plans have been shown to be encoded across the brain [36,41,43]. Second, variance shared across hemispheres may arise from direct anatomical projections [95], which have been shown to be spatially patterned [96]. Interactions across the hemispheres are central to various brain functions such as working memory [48,49], motor control [41,42], sensory perception [40], and spatial attention [56,57].”\n● Reworked discussion paragraph (L399):\n“Underlying the across- and within-area circuit motifs are various cellular and circuit mechanisms, which have been implicated in the shaping of correlations among neurons both across and within brain areas. First, the structure of inhibition is known to shape correlations among neurons, where the influence of inhibition depends on its spatial scale [26,104] and timing relative to excitation [105-108]. Second, changes at the level of the synapse, including synaptic depression [109] and connection strength [110], can also regulate correlated variability. Third, neuronal correlations are influenced by neuromodulators such as acetylcholine [111-113] and norepinephrine [91], which together have wide-ranging modulatory effects including altering the strength of thalamocortical inputs [114] and regulating global brain state [115]. Finally, the strength of excitatory synaptic connectivity has been linked to the strength of correlation [116]. Anatomically, long-range horizontal connections within a brain area [117-119], across brain areas [96], and across the corpus callosum [95], are all weaker than local recurrent connectivity. However, the strength of the detected correlations need not be directly linked to the strength of the anatomical connection between hemispheres or among distant neurons [120]. For example, weak neuromodulatory inputs could drive strong correlations among neurons by influencing the much stronger local circuit mechanisms that exist in local circuitry. Our analyses using pCCA-FA, like any other assessment based solely on spiking activity, do not explicitly distinguish between these cellular and circuit mechanisms. Instead, we focused on the spiking variability yielded by these mechanisms and dissected the components of spiking variability that are shared across versus within brain areas.”\n7. Related to point 6 above, the authors hypothesize that inter-area correlation is driven by a global neuromodulatory signal (such as NE) potentially linked to arousal. If the signal is global, why would it preferentially strengthen correlations between neurons in different brain areas in comparison to correlations between neurons in the same area? The logic here needs additional clarification.\nWe agree with the reviewer that global neuromodulation such as NE is likely to boost statistical interactions among all neurons, both within and across areas. In the pCCA-FA model, the across-area component is defined as variance that is shared with one or more neurons in the other area (regardless of whether that variance is also shared with neurons in the same area). By contrast, the within-area component is defined as variance that is shared exclusively with neurons in the same area. Thus, if global neuromodulation boosts interactions among all neurons, it would show up in the across-area component. We note that this logic was not made clear and appreciate the suggestion to clarify this point. We have added a statement in the text surrounding Fig. 6.\nChanges:\n● L328 in Results for Fig. 6:\n○ Addition: “Recall that the across-area component is defined as variance that is shared with one or more neurons in the other area (regardless of whether that variance is also shared with neurons in the same area).”\n8. Also related to point 6 above, it is thought that more synapses onto cortical pyramidal neurons derive from neighboring neurons in the same brain area than from neurons in distant areas. If synapses from recurrent local circuits do in fact provide the bulk of synaptic input, why then do inputs from remote brain areas drive stronger correlations? It could be useful to consider the correlation results in the context of what is known about the architecture of local versus long-tract connections, and to reconcile apparent discrepancies between the relative strength of connections and the relative strength of correlations.\nWe agree with the reviewer in that the level of interactions we observed across the hemispheres of PFC is seemingly at odds with the strength of anatomical connections. However, we note that the strength of the detected interactions need not be directly linked to the strength of the anatomical connection between hemispheres or among distant neurons (Goulas, et al., 2015). For example, weak neuromodulatory inputs could drive strong correlations among neurons by influencing the much stronger local circuit mechanisms that exist in local circuitry. As noted above in response to point 6, we have added a paragraph to the Discussion that addresses this point.\nMinor\n1. The Methods describe training a neural network to distinguish neural action potential waveforms from noise. Were action potential waveforms from different neurons discriminated? If not, it would be important to state that the analysis was conducted on multi-neuron firing rates.\nThanks for pointing this out. Waveforms from different neurons were not discriminated. Rather, we analyzed multi-unit activity from each electrode on the Utah array. We have now made this clear in the Methods.\nChanges:\n● L536 in Methods section “Preprocessing of neural activity”:\n○ Addition: “The spikes from each electrode on the Utah array were treated as multi-unit activity, and no further differentiation of single units was performed. Thus, each electrode was treated as one neural unit. ”\n2. The paper does an excellent job of introducing the statistical concepts behind the technique. Any additional walk-throughs of the main concepts that the authors can provide would be to the good, to clarify the relationship between the loading matrices, latent variables, and population activity vectors, and how they combine to generate the co-fluctuation in firing patterns illustrated schematically in Fig. 3b.\nWe thank the reviewer for their positive comment. To provide a better walkthrough of the schematics in Figure 3, we have added an additional panel (b). This panel illustrates the concept of a co-fluctuation pattern as a vector of weights (i.e., a column of the loading matrix), providing a link between the parameters estimated by pCCA-FA (loading matrices and latent variables) and the neural activity.\nChanges:\n● Fig. 3b, pasted in context with Fig. 3a and 3c:\n● Caption:\n“b, The activity of each neuron is approximated as a weighted combination of the across-area latent variable (pink) and the within-area latent variable (light and dark blue for areas 1 and 2, respectively). The across-area latent variable (z, same for both areas) is mapped to the activity of each neuron via a pair of across-area co-fluctuation patterns, one for each area (top and bottom, pink). A co-fluctuation pattern is a vector that comprises one weight for each neuron. Each within-area latent variable (z_1 and z_2 for areas 1 and 2, respectively) is mapped to the activity of each neuron in the corresponding area via a within-area co-fluctuation pattern (area 1: top, light blue; area 2: bottom, dark blue). Note that the within-area co-fluctuation patterns map different latent variables (z_1 and z_2) and, thus, are not paired. When there is more than one latent variable of a given type, we simply append analogous terms to the weighted sum.”\n● L197 in Results for Figure 3:\n○ Addition: “pCCA-FA approximates neural activity as a weighted sum of the across-area latent variables (Fig. 3b, pink) and the within-area latent variables (Fig. 3b, light and dark blue for areas 1 and 2, respectively). The co-fluctuation patterns comprise the weights, which are estimated from the neural activity and fixed across all trials. The values of the latent variables can vary from trial to trial.”\nResponse to reviewers\nInteractions across hemispheres in prefrontal cortex reflect global cognitive processing\nMegan E. McDonnell*, Akash Umakantha*, Ryan C. Williamson*, Matthew A. Smith†, & Byron M. Yu†\n* and † indicate equal contribution\nWe would like to thank the reviewers for their time and their helpful comments throughout the review process. We have added one additional citation as requested by Reviewer 1. Below, we have provided a brief point-by-point response (in pink) to each reviewer comment (in black).\nReviewer #1 (Remarks to the Author):\nThe authors have done an excellent job at addressing all the comments. I do not have any further points to add except that a recent Nature paper (Raut at al., 2025) reached a very similar conclusion. It would be useful to discuss how the results of both papers relate, and how pCCA-FA confers any benefit (if any) compared to the time delay embedding used by Rout and colleagues.\nRaut, R. V. et al. Arousal as a universal embedding for spatiotemporal brain dynamics. Nature 647, 454–461 (2025).\nWe thank the reviewer for suggesting this paper. We have added it as a citation in the Discussion section about brainwide arousal, as well as our discussion of targeted dimensionality reduction techniques. The technique applied by Raut et al. specifically aimed to uncover arousal dynamics and relate them to neural activity across the brain (i.e., a supervised method). By contrast, pCCA-FA is an unsupervised method that can be applied to any two simultaneously recorded populations of neural activity for the purpose of identifying variance shared across and within brain areas. Each type of shared variance can then be related to arousal. We agree that both approaches are valid and interesting techniques to uncover latent structure in neural activity, although the goals of each method differ.\nReviewer #3 (Remarks to the Author):\nI would like to thank the authors for the additional work they invested in revising their manuscript. They have answered all my comments to a more-than-satisfactory degree. I am happy to recommend publication of the manuscript in its current form.\nWe thank the reviewer for their remarks, and for their efforts in reviewing this work.\nReviewer #4 (Remarks to the Author):\nThe authors have fully and thoroughly addressed all my prior concerns. They have promoted data to the main text providing strong evidence validating their novel analytical technique. They added control analyses to demonstrate detected correlations depend on the simultaneity of the neural activity. They added control analyses shuffling neurons between hemispheres to show strong impact on detected within-area correlations (thanks to the author for the clarification of the predicted results of this analysis). They provide control analyses showing that differences in the strength of detected correlations within and between areas is not due to the number of free parameters in the models fit to neural activity in the two cases. They have addressed the underlying circuit motifs and synaptic mechanisms that might drive the detected correlations.\nThe revision was remarkably thorough. Each point was considered in depth, with new data provided to address the point.\nThis is a strong and exciting study.\nI look forward to wide adoption of the techniques the authors have developed to detect interactions within and between cortical areas in large neural datasets. The underlying biological questions are important and perhaps understudied in experimental neuroscience – namely how do interactions between neurons that impose correlation on their activity drive computation in distributed systems.\nWe thank the reviewer for their kind words, and for their efforts in reviewing this work."}]

Decision letter