Using a virtual-reality system, we trained three monkeys (monkeys H, N, and S) to reach for a visual target with their nonvisible (proprioceptive) arm while viewing a virtual arm moving in synchrony with a preset spatial visual-proprioceptive (VP) disparity (Figure 1A). On each trial of the experiment, the monkeys were required to initiate the trial by placing their hand on the starting position (blue dot) for 1 s and were instructed not to move. After the initiation period, the starting point disappeared, and the visual virtual arm was rotated; this mismatch arm was maintained for 0.5 s as the preparation period. The reaching target was presented as a ‘go’ signal, and monkeys had to reach toward the visual target within 2.5 s and place their hand in the target area for 0.5 s, referred to as the target-holding period, to receive a reward. Any arm movement during the target-holding period automatically terminated the trial. The proprioceptive drift due to the disparity between visual and proprioceptive inputs was measured at the endpoint of the reach and was defined as the angle difference between the proprioceptive arm and the visual target (the estimated arm) (Figure 1B, see details of animal training and reward in Materials and methods). In addition to this VP conflict (VPC) task, two control experiments were conducted: (i) where the visual and proprioceptive signals were perfectly aligned (VP task) and (ii) where there was only a proprioceptive signal (P task). The procedures of the three tasks (VPC, VP, and P) were identical, except that the visual or proprioceptive information presented to monkeys varied according to the context of the experiment (see Materials and methods). Using a block design, the order of three different blocks (tasks) in each training or recording session was randomized.

Figure 1 with 2 supplements see all

Download asset Open asset

Figure 1—source code 1

Related to Figure 1G.

https://cdn.elifesciences.org/articles/76145/elife-76145-fig1-code1-v2.zip

Download elife-76145-fig1-code1-v2.zip

Figure 1—source data 1

Related to Figure 1G.

https://cdn.elifesciences.org/articles/76145/elife-76145-fig1-data1-v2.zip

Download elife-76145-fig1-data1-v2.zip

We used this paradigm to test the hypothesis of the causal inference process, which predicts how the brain infers and updates hidden structures on the basis of multiple sensory inputs.

First, the BCI model encodes probability distributions over the sensory (visual and proprioceptive) signals and incorporates rules that govern how a prior belief about the sensory causal structure is combined with incoming information to judge the event probability in proprioception (Figure 1C–F). Thus, the monkey’s behavior output (the proprioceptive drift distribution under each disparity) should show the dynamics of integration and segregation, which is the hallmark of causal inference. That is, the drift should increase for small disparities and decrease when the disparities become larger (Fang et al., 2019; Körding et al., 2007).

Second, according to the model, the prior of common source in the current trial should be modulated by the experience of the environmental structure. Thus, at a single trial level, the prior of a common source in the current trial should be updated based on the posterior of a common source from previous trials (Figure 1C and D).

Third, the sensory uncertainty is also proposed to update to maintain consistency with the prior beliefs of the causal structure of the world (French and DeAngelis, 2020). Therefore, the sensory uncertainty should increase when there is a conflict between the proprioceptive and visual signals (e.g., the VPC task) (Figure 1D and E).

To test these hypotheses, we adopted the Bayesian causal inference (BCI) model to assess monkeys’ behavior and investigated whether the neural activities in multiple brain regions correlate to proposed components in the behavior.

To examine whether the monkeys inferred the causal structure during multisensory processing, we first examined the proprioceptive drift as a function of disparity in the VPC task. Overall, the three monkeys showed a very consistent behavioral pattern, with the proprioceptive drift increasing for small levels of disparity and plateauing or even decreasing when the disparity became larger (e.g., exceeded 20°) (Figure 1G; for data on individual monkeys, see Figure 1—figure supplement 1). The BCI model qualitatively explains the nonlinear dependence of drift as a function of disparity. For small disparities, there is a high probability that the proprioceptive and visual signals came from the same source. Hence, the visual information is fully integrated with the proprioceptive information. For large disparities, however, the proprioceptive and visual signals are likely from different sources, leading to a breakdown of integration and consideration of only the proprioceptive information (segregation). In this case, visual information has a weaker weight for perception. Consequently, the effect of disparity on the drift is reduced by shifting integration to segregation. The BCI model quantified the nonlinear dependence between disparity and proprioceptive drift to measure the posterior probability of a common source (P_com), the consequence of causal inference. We fitted the behavioral data using the BCI model. The results showed two signatures of the P_com pattern: (i) the averaged P_com decreased as the disparity increased (Figure 2A, left) and (ii) within each disparity, especially the large ones, the P_com decreased as the proprioceptive drift decreased (Figure 2A, right) (see individual monkeys’ behavior in Figure 1—figure supplement 1).

Figure 2 with 2 supplements see all

Download asset Open asset

Figure 2—source code 1

Related to Figure 2A-D.

https://cdn.elifesciences.org/articles/76145/elife-76145-fig2-code1-v2.zip

Download elife-76145-fig2-code1-v2.zip

Figure 2—source data 1

Related to Figure 2A–D.

https://cdn.elifesciences.org/articles/76145/elife-76145-fig2-data1-v2.zip

Download elife-76145-fig2-data1-v2.zip

More importantly, the model posits that not only the inference of the causal structure is based on visual and proprioceptive inputs but also the subsequent updating of (i) the prior belief of causal structure based on the experience (e.g., probability of a common source in the previous trials) and (ii) the uncertainty of sensory signals for the visual and proprioceptive recalibration (Figure 1C–F). To test these hypotheses, we first implemented the Markov analysis of the prior belief and P_com (see Materials and methods) to see whether the prior probability of a common source (P_prior) in the current trial depended on the previous P_com (Figure 2B). The Markov model included the transition probability of P_prior between the current (nth) and previous (n^th− 1) trial to account for the trial-by-trial variability in spatial drifts observed in the three monkeys (Figure 2B, left). The fit to the model demonstrated that the P_com observed in the nth trial was significantly affected by that in the previous (n^th− 1) trial (Figure 2B, right, Wilcoxon signed-rank test, p<0.001), indicating that the P_com was computed based on both P_prior from the previous trial and the sensory inputs, with their disparity, from the current trial. Note that the transition probabilities (P_(C=1|C=1) and P_(C=1|C=2)) remained relatively high (larger than 0.8 in three monkeys) because overall, the number of high P_com trials was much more than low P_com trials in either training or recording sessions. This was consistent with high baseline P_prior in three monkeys (Supplementary file 1).

We next examined whether the sensory representation is updated to maintain consistency with the causal structure of the environment. That is, the estimates of physical arm locations should tradeoff systematically depending on the current common-source belief (e.g., P_com in different tasks: VP, P, and VPC). For example, when the monkey incorrectly infers that the visual and proprioceptive arms come from the same source when a disparity is presented, the uncertainty of proprioception should increase to ‘explain away’ the conflict between the two inputs. According to this idea, since the block design in the current experiment resulted in P trials (in the P task) sometimes following the VPC task and other times following the VP task, we then reasoned that because the overall P_com was lower in the VPC task than in the VP task, the uncertainty of proprioception (i.e., the distribution of proprioceptive drifts in the P trials) would be larger after the VPC task than after the VP task. We analyzed the drift variation in P trials and found that, in the early trials (first third of each P block), the uncertainty of P trials following the VPC task was significantly larger than that following the VP task (Figure 2C, right, Wilcoxon signed-rank test, p=0.012). The increase in the uncertainty of proprioception was recovered in the late trials (last third of each P block), evident by a significant difference in the uncertainty between early and late P trials (Figure 2C, right, Wilcoxon signed-rank test, p=0.035). The decrease in the uncertainty of proprioception was reasonable, as the tradeoff effect in the VPC task gradually recovered.

Furthermore, we hypothesized that if a tradeoff of sensory representation occurs during the process of causal inference, the tradeoff would also affect the uncertainty of VP integration in both VP and VPC tasks. We examined the distribution of proprioceptive drifts using the trials with 0° disparity in the VPC task, in which the V and P information were congruent, and compared it with the distribution in the VP task. As predicted, we found that the variance of the proprioceptive drift was significantly larger in the VPC task than in the VP task (Figure 2D, right, Wilcoxon signed-rank test, p<0.001). Note that the difference between VP and VPC (0°) tasks could not be explained by the divergence of eye fixation positions (Figure 2—figure supplement 2). As a control, we also investigated whether the mean of drift, representing the perceptual accuracy of the proprioceptive arm, was affected by the causal structure of the environment. We found there was no significant difference between the mean of drift for P trials following the VPC task and that following the VP task in both early parts (Figure 2—figure supplement 1, left, Wilcoxon signed-rank test, p=0.37, false-discovery rate [FDR] corrected) and late parts (Figure 2—figure supplement 1, right, Wilcoxon signed-rank test, p=0.37, FDR corrected). Besides these, we also found that the mean of proprioceptive drift was not updated in the VPC task compared with the VP task (Figure 2—figure supplement 1, right, Wilcoxon signed-rank test, p=0.29). Thus, these results supported the notion of a tradeoff in proprioception according to causal inference environments; that is, sensory representation’s uncertainty, not accuracy, is updated dynamically based on the task environment (P_com).

To summarize the above-described behavioral results, we found that monkeys’ proprioceptive drift shows a nonlinear dependency on the disparity between proprioceptive and visual input, which was well explained by the causal inference model. Second, we showed that the P_com integrated with VP sensory inputs and is updated by previous experience on a trial-by-trial basis. Third, to maintain a consistency of causal inference, sensory uncertainty, reflected by the variance of proprioceptive drift, is updated in the inference along with the change of P_com. Taken together, we established the behavioral paradigm in which monkeys infer the hidden cause by integrating prior information and sensory inputs while dynamically updating both P_com and sensory representation. The behavioral responses of the monkeys enabled us to examine the underlying neural mechanisms and functional circuits.

Previous studies showed that the premotor and parietal cortices were highly involved in body representation and multisensory perception (see reviews in Blanke, 2012; Graziano and Botvinick, 2002). In monkeys, bimodal neurons with visual and somatosensory receptive fields were found in both premotor (including F2vr in dorsal premotor and F4/F5 in ventral premotor) and posterior parietal cortices (including area 5 and area 7) (Fogassi et al., 1999; Graziano et al., 2000; Graziano and Gross, 1993; Graziano and Gross, 1998; Graziano et al., 1994). Specifically, ventral premotor neurons responded to visual stimuli in the space adjacent to the arm (Graziano and Gross, 1998; Graziano et al., 1994). The bimodal neurons in the parietal cortex (area 5 and area 7) showed to respond to both the real arm position and the seen position of a dummy arm (Graziano et al., 2000), which have a significant projection of the premotor cortex (Graziano and Gross, 1998). Consistently, human fMRI studies found that the posterior parietal and premotor (dorsal and ventral) cortices selectively respond to visual stimulation near the hand (Brozzoli et al., 2011) or the dummy hand near one’s corresponding hand (Blanke et al., 2015; Ehrsson et al., 2004). A human MEG study also revealed that the activities in the prefrontal and intraparietal sulcus were related to the causal inference computation in visual-auditory integration (Cao et al., 2019; Rohe et al., 2019). Therefore, we determined to record from two brain regions, the premotor cortex (dorsal and ventral, 412 neurons) and parietal cortex (area 5 and area 7; 238 neurons), in the three monkeys performing the reaching tasks (Figure 3A, for details, see Materials and methods). We first examined whether neurons in the premotor and parietal cortices during the target-holding period (Figure 3B) were selective to basic task components, including condition (VP or P), arm location, and visual disparity. In the premotor cortex, 40% (163/412) of neurons were selective to condition, 23% to arm location, and 37% to visual disparity (Figure 3—figure supplement 1A, upper panel). In the parietal cortex, 35% (83/238) of neurons were selective to condition, 27% to arm location, and 31% to visual disparity (Figure 3—figure supplement 1A, lower panel, ANOVA, main effect, p<0.05). We also examined the neural representations of the visual and proprioceptive arm locations in each trial during the target-holding period in the VPC, VP, and P tasks, measured by a bias-corrected percent explained variance (ωPEV) (Figure 3C). Both brain regions conveyed vital information about the arm location in the three tasks. In the VP and P tasks with no VP disparities, both premotor and parietal cortices showed similar visual and proprioceptive arm information (Figure 3C). However, when disparities were introduced in the VPC task, the premotor cortex showed a more robust signal for visual arm information (Figure 3C). In contrast, the parietal cortex showed stronger signals for information related to the proprioceptive arm (Figure 3C).

Figure 3 with 3 supplements see all

Download asset Open asset

Figure 3—source code 1

Related to Figure 3C, D, F, G and H.

https://cdn.elifesciences.org/articles/76145/elife-76145-fig3-code1-v2.zip

Download elife-76145-fig3-code1-v2.zip

Figure 3—source data 1

Related to Figure 3C, D, F, G and H.

https://cdn.elifesciences.org/articles/76145/elife-76145-fig3-data1-v2.zip

Download elife-76145-fig3-data1-v2.zip

Next, to define causal inference response in the VPC task at the single-neuron and single-trial levels, we utilized the VP and P tasks to characterize neural responses, as these tasks involve expected stereotypical behaviors in the two extreme regimes: full integration and segregation. Thus, neurons that are more active during the P task are likely candidates for ‘segregation (P) neurons’, which exhibited increased activity under the large disparities in the VPC task (Figure 3D). By contrast, neurons that are more active during the VP task reflect a preference for integrating congruent VP information and, hence, constitute a natural candidate for ‘integration (VP) neurons’ (example in Figure 3—figure supplement 2). We then implemented a linear probabilistic model which combined how the neural response pattern aligned with the VP and P response profiles and used this model to implement a probabilistic decoding analysis to calculate the probability of VP or P (VP weight = P_vp/[P_vp + P_p]) based on the firing rate in each trial (Figure 3E; also see Materials and methods). Thus, a larger VP weight for a single trial denotes a higher probability of integration (high P_com). We first focused on the target-holding period in a trial, as the neurons could well display their spatial tunings when monkeys holding their arms on the target. We found that both premotor and parietal cortices carry information about P_com at the single-neuron (Figure 3F; the same example neurons in Figure 3D) and population levels (Figure 3G; see Materials and methods) during the target-holding period. That is, the VP weight of the neuron or population progressively decreased along with the disparity, and in trials with large disparity (e.g., 35° and 45°), the neuron(s) had a higher VP weight when the drift was large (i.e., the monkey integrated the visual information; thus, a high P_com predicted by the BCI model) and shifted gradually toward higher P weights when the drift shifted to 0 (i.e., the monkey segregated the visual information; thus, a low P_com predicted by the BCI model). The VP weight was highly correlated with the P_com from behavior (Figure 3H). Note that the premotor cortex had a slightly higher proportion of causal inference neurons (11.7%) than the parietal cortex (7.6%, Pearson chi-square test, χ²=2.33, p=0.063).

As neuronal activities in the premotor and the parietal cortices are reported to correlate with the eye position in the reaching task (Buneo and Andersen, 2006; Pesaran et al., 2006), one might ask whether the P_com signals can be explained by the eye position. However, the result showed that the VP weights in the population could not be predicted by eye fixation positions during the target-holding period (Figure 3—figure supplement 3).

We next focused on the overall populations of neurons in both regions and asked whether and how their population states reflect the uncertainty of causal structure, P_com. We were guided by the results from single-neuron analyses during the target-holding period described above, in which neurons responsive to high P_com (prefer integration) are more likely to show neural tuning similar to that during the VP task, and neurons responsive to low P_com (prefer segregation) show a tuning profile similar to that in the P task. We thus hypothesized that neural components or subspaces embedded in the population activity represent the dynamic change in the coding of P_com in the VPC task, which would lie between the components representing the VP and P profiles. Furthermore, the computation of P_com in the BCI model is determined by the relation and disparities between the visual information from the artificial arm and proprioceptive information from the monkey’s actual arm. In other words, according to the model, the causal inference can be constructed before the visual target appears, and the participant uses this information to guide the reach. We thus further hypothesized that the dynamics of the population states also reflect the P_com during the preparation period, during which there is no motor planning or preparation.

Thus, we grouped trials from each neuron into high and low P_com classes according to the drift under each disparity (high, the top third of the trials [in red]; low, bottom third of the trials [in blue]) (Figure 4A). We conducted demixed principal component analysis (dPCA) to visualize any neural component that represents the P_com in the VPC task in relation to that in the VP and P tasks (see Materials and methods). dPCA decomposes population activity into a set of dimensions that each explain the variance of one factor of the data (Kobak et al., 2016). We included the factors of time, arm location, and P_com (Figure 4B). In the analysis, VP and P trials were included, which served as the templates of integration and segregation, respectively. As shown in the schema (Figure 4B), if the decomposed neural components indeed represent the P_com, the population activity of high and low classes in this subspace should lie between that of the VP and P classes and the four classes (high, low, VP, and P) should be separated from each other. The dPCA results indicated that the P_com components, unrelated to the arm location, represented 29.9% and 20.5% of the total firing rate variance in the premotor and parietal cortices, respectively (Figure 4C, in red). Notably, the activity in P_com dimensions seems consistent with our hypothesis, demonstrating the dynamics of P_com between integration (VP) and segregation (P). In addition, compared to the activity in the parietal cortex, the neural trajectories of the premotor populations showed an earlier divergence in P_com dimensions (Figure 4D).

Figure 4 with 3 supplements see all

Download asset Open asset

Figure 4—source code 1

Related to Figure 4C–F.

https://cdn.elifesciences.org/articles/76145/elife-76145-fig4-code1-v2.zip

Download elife-76145-fig4-code1-v2.zip

Figure 4—source data 1

Related to Figure 4C–F.

https://cdn.elifesciences.org/articles/76145/elife-76145-fig4-data1-v2.zip

Download elife-76145-fig4-data1-v2.zip

To further quantify their dynamics statistically, we trained a linear support vector machine (SVM) using pooled activities in each brain region throughout the entire trial. The dynamic decoding results showed that the P_com information is correctly predicted by neuronal population activities in both areas after target onset but is decoded only by premotor neurons during the preparation period when there was no visual target or motor preparation (Figure 4E, cluster-based permutation test, p<0.05). Randomization test confirmed the time difference that the P_com information occurred significantly earlier in the premotor cortex than the parietal cortex (Figure 4—figure supplement 3, randomization test, p<0.01, see Materials and methods). This may suggest that the premotor cortex is where causal inference is computed and sends the information to the parietal cortex during the reaching period.

Next, we tested the relationship between the population activities in the two areas. We performed a joint peri-event canonical correlation (jPECC) analysis, which detects correlations in a ‘communication subspace’ between two brain regions (Steinmetz et al., 2019). In brief, we conducted a canonical correlation analysis for every pair of time points containing the population neural firing rates from the two regions. If the shared neural activity emerges at different times in the two areas, that is, activity in one region potentially leads to activity in the other, then we should observe a temporal offset between them. The jPECC results revealed a significant time lag for activity correlations between premotor and parietal areas in P_com dimensions (Figure 4F, cluster-based permutation test, p<0.05), suggesting a potential feedback signal of P_com from the premotor cortex to the parietal cortex. As a control, we performed the same procedure with misalignment trials (see Materials and methods) to exclude the probability that the observed time lag resulted from the intrinsic temporal property of neuronal activities in these regions. There was no significant time lag between premotor and parietal areas when the trials were misaligned (Figure 4—figure supplement 1).

The behavioral experiments showed that the P_com could be updated by previous sensory experience on a trial-by-trial basis. To test the effect of the previous P_com on the causal inference in each trial, we examined neural activities during the baseline period in the VPC task before a disparity in the visual and proprioceptive arm was introduced (Figure 5A). We again classified the trials according to high and low P_com. Figure 5A depicts the results from an example premotor neuron, showing that during the baseline period, the neural activity exhibited selectivity toward the previous trial’s P_com, and at the same time, its neural trajectories in high and low prior classes lay between the VP and P templates. Of 412 neurons in the premotor cortex, 29 (7.0%) showed such selectivity in the previous trial (Figure 5—figure supplement 1).

Figure 5 with 2 supplements see all

Download asset Open asset

Figure 5—source code 1

Related to Figure 5A–F.

https://cdn.elifesciences.org/articles/76145/elife-76145-fig5-code1-v2.zip

Download elife-76145-fig5-code1-v2.zip

Figure 5—source data 1

Related to Figure 5A–F.

https://cdn.elifesciences.org/articles/76145/elife-76145-fig5-data1-v2.zip

Download elife-76145-fig5-data1-v2.zip

To further test the relation between baseline neural activity and behavior quantitatively, we examined whether the population activities of these neurons can predict the P_com from previous trials. We trained an SVM using pooled activities across recording sessions. The previous P_com was only correctly decoded from the baseline activity in the premotor cortex (Figure 5B, cluster-based permutation test, p<0.05). Moreover, only recent experience (nth−1 trial) had a significant impact on the current trial (Figure 5C, permutation test, p<0.001).

As both P_prior and P_com were represented in premotor neural activities, we wanted to examine their relationship in the neural states. We first found that very few neurons responded to both information types (see Figure 5—figure supplement 1). We then hypothesized that P_prior and P_com might be represented independently at a population level. To validate this hypothesis, we conducted PCA on the population activities during baseline and target-holding periods for P_prior and P_com, respectively. If they are independent, the subspaces of P_prior and P_com will be near orthogonal, and the PCs of P_prior and P_com will capture little variance from each other (Elsayed et al., 2016). To quantify this, we projected the P_prior data onto the P_com subspace to calculate the percent variance explained by the P_com PCs and repeated the same procedure for the P_com data (Figure 5D). The results show that the top 10 P_prior PCs captured very little P_com variance; similarly, the top 10 P_com PCs captured very little P_prior variance (Figure 5E). These results support the hypothesis that the two information types are represented independently in the premotor cortex. However, such independence between P_com and P_prior could also be caused by their different temporal structures in the task. Thus, we examined their neural dynamics within a trial. Figure 5F shows the time course of decoding results of prior and posterior information, where the P_prior quickly decreased after the disparity onset. At the same time, the P_com information increased and was retained until the end of the trial. These results demonstrated the dynamics in the computation of causal inference, where the information from the last trial is only preserved transiently and then used to integrate with sensory inputs to generate P_com information.

Finally, we investigated the neural activities associated with updating sensory uncertainty. The behavior results revealed a significantly greater uncertainty of proprioception in VP trials in the VPC task (low belief of a common source) than in the VP task (high belief of a common source) (Figure 2D). We hypothesized that the sensory signals, which were used to make causal inference, in turn, updated their neuronal tunings to match inferred causal structure. We first examined the difference in neural tuning for arm location using the VP trials in the VP and VPC (VPC (0°), trials with no disparity) tasks. To test whether the tuning functions of arm location selective neurons changed between the VP condition and VPC (0°) condition at the single-neuron level, we fitted the tuning curve with the von Mises distribution by using the neuron response in different arm locations (five levels: [−30°, −20°, 0°, 20°, and 30°]) for these two conditions respectively (see Materials and methods). We found that the averaged firing rates during the target-holding period under the VP condition were higher than that under the VPC condition (0°) in the parietal cortex (Figure 6—figure supplement 1, left, Wilcoxon signed-rank test, p=0.017) but not in the premotor cortex (p=0.71). The gain index under VP condition were higher than the VPC condition (0°) in the parietal cortex (Figure 6—figure supplement 1, middle, Wilcoxon signed-rank test, p=0.0016, FDR corrected) but not the premotor cortex (p=0.11, FDR corrected). Figure 6A (right) shows an example neuron from the parietal cortex tuned to the center (0°) of arm location in the VP task, and the tuning range/uncertainty of the arm location was broader/lower in the VPC task. Here, for visualization purposes, we selected the time point when this neuron demonstrated the highest difference of ωPEV in the VP trials between VP and VPC tasks for the tuning calculation (Figure 6A, left, peak delta ωPEV). The averaged dynamic spatial selectivity of all neurons revealed a significant decrease of the total spike rate variance explained by the arm location in the parietal cortex but not in the premotor cortex (Figure 6B, cluster-based permutation test, p<0.05). Note that the updating of sensory uncertainty was not correlated with the uncertainty of eye position between VP and VPC (0°) tasks (Figure 2—figure supplement 2).

Figure 6 with 1 supplement see all

Download asset Open asset

Figure 6—source code 1

Related to Figure 6A–C.

https://cdn.elifesciences.org/articles/76145/elife-76145-fig6-code1-v2.zip

Download elife-76145-fig6-code1-v2.zip

Figure 6—source data 1

Related to Figure 6A–C.

https://cdn.elifesciences.org/articles/76145/elife-76145-fig6-data1-v2.zip

Download elife-76145-fig6-data1-v2.zip

Furthermore, at the population level, we performed the SVM decoding analysis of arm locations and found that only the parietal cortex showed a significantly decreased decoding accuracy in the VPC task (Figure 6C, cluster-based permutation test, p<0.05). We also confirmed that the change of decoding accuracy in the parietal cortex was significantly larger than the change in the premotor cortex (two-way ANOVA, Condition (VP and VPC (0°))×Region (premotor and parietal), significant interaction effect, p<0.05).

All animal procedures were approved by the Animal Care Committee of the Center for Excellence in Brain Science and Intelligence Technology, Institute of Neuroscience, Chinese Academy of Sciences (Permit Number: CEBSIT-2020034), and were described previously in detail (Fang et al., 2019). Three male adult rhesus monkeys (Macaca mulatta; Monkeys H, N, and S, weighting 6–10 kg) participated in the experiment. During the experiment, the monkeys were seated comfortably in the monkey chairs, and their heads were fixed. All monkeys were implanted with chambers for recordings.

Some of the following methods are similar to those previously published (Fang et al., 2019).

The monkeys were seated in front of a chest-height table on which a lab-made virtual reality system was placed (Fang et al., 2019). During the experiment, the monkey’s left arm (and the right arm in the case of Monkey H, who was right-handed) was placed in the system and blocked from sight. A CCD camera (MV-VEM120SC; Microvision Co., China) captured the image of the monkey’s arm reflected in a 45° mirror. This image was projected to the rear screen by a high-resolution projector (BenQ MX602, China). Therefore, when the monkey looked in the horizontal mirror suspended between the screen and the table, the visual arm image appeared to be its real arm on the table. The lower edge of the screen was aligned to the table edge. The monkey’s trunk was close to the edge of the table, and the left shoulder was aligned with the midline of the screen. Using the OpenCV graphics libraries in C++ (Visual Studio 2010; Microsoft Co., Redmond, WA, USA), the arms image and the visual target were generated and manipulated. Using CinePlex Behavioral Research Systems (Plexon Inc, Dallas, TX, USA), sampled at 80 Hz, the hand position was tracked and recorded. The tracking color marker was painted onto the monkey’s first segment of the middle finger, which was not visible after adjusting the light exposure settings of the video.

The monkey was trained to report its proprioceptive arm location by reaching for a target in a VPC causal inference task (Figure 1A; Fang et al., 2019). The monkey initiated a trial by placing its hand on the starting point (a blue dot with a 1.5 cm diameter) for 1000 ms and was instructed not to move. After the initiation period, the starting point disappeared, and the visual arm was rotated (within one video frame, 16.7 ms) for the VPC task. The rotation was maintained for 500 ms (the preparation period). After that, the reaching target was presented as a ‘go’ signal. The monkey had to reach the target (chosen from T1 to T5 randomly trial by trial [Figure 1A]) within 2500 ms and hold its hand in the target area (see as follows) for 500 ms to receive a drop of juice as the reward. Any arm movement during the target-holding period automatically terminated the trial. The rotated arm was maintained throughout the entire trial along with the arm movement. The intertrial interval (ITI) was ~1.5–2 s, after which the monkey was allowed to start the subsequent trial. During the ITI, the visual scene was blank. Under the VPC task, across trials, the visual arm was randomly presented with a disparity of 0°, ±10°,± 20°, ±35°, or ±45° (+, clockwise [CW]; −, counterclockwise [CCW] direction) from the subject’s proprioceptive arm, with its shoulder as the center point. The starting point was fixed 25 cm away from the monkey’s shoulder. The target position was selected randomly trial by trial from one of five possible positions located on an arc (a ±4° jitter was added to the original position trial by trial to ensure the monkey did not perform the task by memorizing all the target positions).

Besides the VPC task, the monkey was also instructed to perform a VP congruent and P task during the recording session. The only difference between the VPC and VP task was that during the entire trial under the VP task, the visual arm was always congruent with the proprioceptive arm. The only difference between VP and P tasks was that during the single trial for the P task, the visual arm information was blocked starting from the onset of the preparation period.

Each VPC block contained 55 trials in which the nine disparities and five targets were randomly combined. Each VP and P block had 27 trials in which five targets randomly occurred in every single trial. In one recording session, typically, one or two P blocks were given first to ensure that the monkey performed the task with its proprioceptive arm, and then in the following blocks, VP, P, and VPC tasks were randomly mixed. One recording session contained more than three VP and P blocks and more than eight VPC blocks.

To ensure the monkeys indeed performed the reaching-to-target task with their proprioceptive hand, under the VPC task, the reaching target area (with reward) was defined as follows: the radial distance from the hand to the center of the target was less than 5 cm to ensure that the monkey did reach out to the target; with the target as the center, the azimuth range was set from [−7 (8 for some sessions, same below) + rotation degree/disparity] to +7° when the rotation degree was negative (counter-clockwise), and from –7° to [+7 + rotation degree/disparity] when the rotation degree was positive (clockwise). As shown in Figure 1B (green zone), the reward area ensured the monkey performed the task rationally and without visual feedback. That is monkey’s reaching position between two extreme conditions: one is that the monkey reaches the target purely relying on the visible arm (the drift is equal to the disparity); the other is that the monkey relies on the proprioceptive arm (the drift is equal to zero). Only the correct trials (when the monkey’s arm was located within the reward zone) were used in the subsequent analysis.

Extracellular single-unit recordings were performed as described previously (Fang et al., 2019; Merchant et al., 2013) from three hemispheres in three monkeys. Briefly, under strictly sterile tasks and general anesthesia with isoflurane, a cylindrical recording chamber (Crist Instrument Co., Inc, Hagerstown, MD, USA) of 22 mm diameter was implanted in the premotor cortex and the parietal cortex (area 5 and area 7). We collected the structural magnetic resonance images (MRI) of three monkeys (3T, Center for Excellence in Brain Science and Intelligence Technology, Institute of Neuroscience, Chinese Academy of Sciences), while they were in an MRI-compatible Horsley-Clarke stereotaxic apparatus. The location of the recording chamber on each animal was determined by the atlas with the origin at the Ear Bar Zero (Saleem and Logothetis, 2012). The centers of implanting recording chambers were [right: 20.0 mm; forward: 10.0 mm] for the premotor cortex in Monkey N, [left: 21.9 mm; forward: 24.9 mm] for the premotor cortex in Monkey H, [right: 14.7 mm; forward: 1.1 mm] for the parietal cortex in Monkey N, and, [left: 17.0 mm; forward: 3.5 mm] for the parietal cortex in Monkey S. During the recording session, glass-coated tungsten electrodes (1–2 MΩ; Alpha Omega, Israel) were inserted into the cortex via a guide tube using a multi-electrode driver (NAN electrode system; Plexon Inc, Dallas, TX, USA). All isolated neurons were recorded regardless of their activity during the task, with the recording locations varying from session to session. At each location, the raw extracellular membrane potential was sampled at 40 kHz. On-line raw neural signals were processed offline to obtain a single unit by Offline Sorter (Plexon Inc, Dallas, TX, USA). All spike data were re-sorted using off-line spike sorting clustering algorithms (Offline Sorter, PCA) (Merchant et al., 2013). With manual adjustments, only well-isolated units were considered for further analysis (signal-to-noise is larger than 3). The sorted files were then exported in MATLAB format for further analysis in MATLAB (Mathworks, Natick, MA, USA) and Python (The Python Software Foundation).

All statistical analyses were implemented with scripts written in MATLAB or Python. In the premotor cortex, 412 neurons were recorded from two monkeys (231 neurons from Monkey H and 181 neurons from Monkey N); in the parietal cortex (area 5 and area 7), 238 neurons were recorded from two monkeys (116 neurons from Monkey N and 122 neurons from Monkey S). As all monkeys’ behavior and model fitting results were similar, for all analyses, data were combined across monkeys. All related statistics are reported in the figure legends.

In all cases, we used the Benjamini-Hochberg procedure (Benjamini and Hochberg, 1995) to control FDR at an α=0.05 level, as follows. The p-values of a given set of hypothesis tests were sorted in ascending order, {p₁, p₂, …, p_n}, and we found the first rank $i_{\alpha }$ such that ${\displaystyle p_{i_{\alpha }}\le i_{\alpha }\times 0.05/n}$ . Then we considered tests to be significantly above chance (rejecting null hypotheses) for all ${\displaystyle p<p_{i_{\alpha }}}$ .

To estimate continuous time-dependent firing rates, timestamps of spiking events were resampled at 1 kHz and converted into binary spikes for single trials. Spike trains were then convolved with a symmetric Hann kernel (MATLAB, MathWorks),

where A is a normalization factor ensuring the sum of the kernel values equals 1. Window width L was set to 300 ms. Single neurons were included in the analysis only if they had been recorded for a full set of tasks (VP, P, and VPC tasks with nine disparities: 0°, ±10°, ±20°, ±35°, and ±45°).

Peri-stimulus time histograms (PSTHs) were then calculated for four periods of interest in a trial: (i) the baseline period (500 ms before the onset of visual arm rotation), (ii) the preparation period (500 ms after the onset of the visual arm rotation), (iii) the target-onset period (1000 ms after the onset of target onset), and (iv) the target-holding period (500 ms after the onset of target-holding). To smooth the firing rate at each time point, the neural firing rate was calculated by averaging in sliding windows (window size, 400 ms; step size, 100 ms) in a single trial (Fried et al., 2011; Gu et al., 2016), resulting in 22 time bins of mean firing rate for every single trial for subsequent dynamic analysis.

To examine whether neurons in the premotor and parietal cortices during the target-holding period were selective to basic task components, including condition (VP or P task), arm location, and visual disparity. For each neuron, we conducted a two-way ANOVA in two datasets. One dataset contains the VP and P tasks (condition (two levels: VP and P tasks)×arm location (five levels: [–30°, –20°, 0°, 20°, and 30°]); the response variable is the mean firing rate during the holding period of each neuron). If a main effect of condition (or arm location) in the two-way ANOVA was found (p<0.05), this neuron was classified as a condition (or arm location) selective neuron. The other dataset is the VPC task (the visual disparity (nine levels: [–45°, –35°, –20°, –10°, 0°, 10°, 20°, 35°, 45°])×target position (five levels: [–30°, –20°, 0°, 20°, and 30°]); the response variable is the mean firing rate during the holding period of each neuron). If the main effect of visual disparity in the two-way ANOVA was found (p<0.05), this neuron was classified as a visual disparity selective neuron.

To investigate whether the tuning functions of arm location selective neurons (Figure 3—figure supplement 1) changed between VP condition and VPC (0°) condition at single-neuron level, we fitted the tuning curve with a reduced von Mises function by using the neuron response in different arm location (five levels: [−30°, −20°, 0°, 20°, and 30°]) for these two conditions separately. Here, VPC (0°) condition represents the trials in VPC condition where the disparity equals to 0. And the fitting function was defined as:

where b is the spontaneous firing rate of the neuron, a is defined as the gain index, and μ is preferred arm location. $fr\left(x\right)$ represents the firing rate when the arm location is x. We analyzed the spontaneous firing rate of the neuron, gain index, and preferred arm location between VP condition and VPC (0°) condition in both premotor and parietal cortices (Figure 6—figure supplement 1).

To measure the representation of a single neuron for causal inference on a single trial, the probability that a single neuron would integrate or segregate the sensory information on a single trial was calculated (Fang et al., 2019). The basic assumption here is that in a single trial under the VPC task, if the neuron is more inclined to represent integrated information, then its firing rate will be closer to its response under VP tasks and farther away from the response under P tasks, and vice versa. The normalized weight of integration (VP weight) was calculated as follows:

(1) First, obtain the neuron response to the arm position under P and VP tasks and fit the von Mises distribution to get the tuning curve.

(2) Under VPC tasks, obtain the current visual arm and the real arm positions, and at the same time, obtain the neuron’s firing rate when the arm is in the corresponding position under VP and P tasks, ${\lambda }_{VP}$ , and ${\lambda }_P$, respectively.

(3) The VP and P templates can be generated through the Poisson distribution:

(4) According to the corresponding probabilities, ${Pr}_{VP}$ and ${Pr}_P$ in the two templates are obtained, and the integration weights for this neuron in the VPC task can be obtained through standardization:

To quantitatively describe whether a single neuron is encoding causal inference, the correlation between P_com and VP weight is calculated. The logic is as follows: the P_com can be used to measure the degree of integration or segregation of sensory information at the behavioral level, whereas VP weight can measure this characteristic at the electrophysiological level. Therefore, if a neuron is performing causal inference, there should be a significant positive correlation between the P_com and VP weight for the corresponding behavior. Neurons that (i) respond to VP/P tasks and (ii) for which P_com and VP weight are significantly positively correlated in the final holding stage are called causal inference neurons. The specific algorithm was as follows:

1. First, obtain neurons with significant selectivity under VP and P tasks (condition selective neuron, see Materials and methods: Task selective neurons).

2. According to proprioception drift, all trials were divided into 29 classes. Continuous drift values were grouped into nine clusters: < –35°, [–35° –25°], [–25° –15°], [–15° –6°], [–6°+6°], [+6°+15°], [+15°+25°], [+25°+35°], >+35°. To be noticed, ±6° covers approximately 99% of drift distribution under the VP and P task. Thus, for the disparity of 0°, there was only one cluster [–6°+6°]. Since the distribution of drift becomes wider (higher variance) the larger the disparity, the more clusters would be assigned for the big disparity. For example, for the disparity ±45°, there were five clusters of drifts. P_com and VP weight were assigned for each class by averaging all trials within it. The Pearson correlation coefficient was then calculated between P_com and VP weight. If the P_com and VP weight were correlated significantly and positively (p<0.05 and r>0), the neuron was called a causal inference neuron.

To visualize the VP weight pattern at the brain region level, the VP weight of each trial of a single neuron under VPC tasks was calculated and then divided into 29 clusters as described above. Then, the bootstrap method was used to randomly select 50 trials from each cluster for averaging. This was repeated 50 times to obtain the VP weight (50×29) of a neuron for visualization. This results in a 50 × 29 × N matrix, where N indicates the number of neurons in each brain region (all neurons were used). The trial corresponding to each neuron was averaged to obtain a 50×29 matrix. The VP weights of a brain region were visualized in a heatmap.

To characterize the dynamic representation of the P_com in the entire session, all trials in a recording session were divided into high P_com trials and low P_com trials based on the relative proprioception drift (RD). Each trial’s relative proprioception drift (RD = drift/disparity) was calculated. The basic idea was that the larger the P_com, the more likely the monkey would integrate the visual and proprioceptive information, and the corresponding RD is closer to 1. The top third and bottom third of the trials were designated the high P_com class and the low P_com class, respectively. These grouping methods were verified by the dPCA.

The method for dPCA was adopted from that published in a previous study (Kobak et al., 2016). Time, target position/arm location (−30°, −20°, 0°, 20°, and 30°), and P_com (VP, P, high P_com, and low P_com) were combined to obtain the marginalized covariance matrix of the three. The neurons whose trial number was not less than five under a single condition were selected for dPCA. Population activity was then projected on the decoding axes and ordered by their explained total variance for each marginalization.

The percentage of explained variance (PEV) (Buschman et al., 2011) was used to measure the basic task components encoded by a single neuron, in which PEV reflected the degree to which the variance of a single neuron can be explained for a specific task component. Generally, PEV can be expressed as a statistical value of ${\eta }^2$ , that is, the variance ratio between groups to the total variance. As the statistical value of ${\eta }^2$ has a strong positive bias for a small sample, the unbiased ${\omega }^2$ statistical value (ωPEV) (Olejnik and Algina, 2003) was used.

To evaluate the information about the locations of the proprioceptive arm, visual arm, and estimated arm encoded by a single neuron in the VPC task, an analysis of covariance was used to decompose the variance, and the ωPEV was calculated. In detail, for a single neuron, ωPEV was calculated for each type of arm when setting the other two types of arm locations as covariates. The whole reaching space was divided into 11 parts from −45° to 45° to transform it from a continuous variable to a discrete variable. A nonparametric Wilcoxon rank-sum test was used for unpaired data for statistical significance analysis comparing two brain regions.

The ωPEV was calculated in each time bin to characterize the temporal dynamics of neural information under VP and VPC (0°) tasks. The baseline was defined as the period 500 ms before the onset of visual arm rotation. A one-sided, paired Wilcoxon signed-rank with FDR correction determined the time bins significantly different from the baseline. The time bins showing significant differences between VP and VPC (0°) tasks were determined by a cluster-based permutation test (Gramfort et al., 2013).

To test the relationship between population activities in the two brain regions, the jPECC method described in a previous study (Steinmetz et al., 2019) was utilized. First, the neuronal responses in two brain regions under the same behavior conditions, namely, high P_com and low P_com, were aligned. Then, a PCA was conducted across time and trials to reduce the dimensionality to obtain the first 10 PCs for each brain region. The trials were then divided into 10 equal parts (training set and testing set) for cross-validation (10-fold cross-validation). The PCs of the training set of each brain region were used to perform canonical correlation analysis to obtain the first pair of canonical correlation components (L2 regularization, λ=0.5). Then, the PCs of the testing set from each brain region were projected onto the first pair of canonical correlation components, and the correlation was determined by the Pearson correlation coefficient between these projections from each region. This analysis was performed for each pair of time bins to construct a cross-validated correlation coefficient matrix. Fifty trials for each group (high P_com and low P_com) from each brain region were randomly selected by bootstrapping in this analysis. Finally, a heatmap was obtained by averaging the correlation coefficient matrix repeated 1000 times.

To quantify the lead-lag relationship of information exchange between brain regions, an asymmetric index was calculated by diagonally slicing the jPECC matrix from +300 ms to +300 ms relative to each time point (Steinmetz et al., 2019). For time point t, the average correlation coefficient across the left half of this slice (i.e., the average along a vector from [t − 300, t + 300] to [t, t]) was subtracted from the right half of this slice (from [t, t] to [t + 300, t − 300]) to yield the asymmetry index. To test the leading significant time point across brain regions, the data from neurons in these brain regions were exchanged, and the above-described analysis was repeated 1000 times to obtain the null distribution of the asymmetric index. Then, a cluster-based permutation test was performed to test whether the symmetric index was significantly greater than the chance level (Gramfort et al., 2013).

To further exclude the possibility that the observed lead-lag relationship resulted from the intrinsic properties of neuronal activities rather than the encoded information in these regions, all trials in each brain region were shuffled to ensure that the inter-region trials were not aligned. Then, the analysis was repeated as described above to obtain the asymmetric index.

Note that, due to the limitations of the asynchronous recording (the premotor and parietal neurons were grouped from different individual animals and only Monkey N was recorded in both areas), further studies are required to clarify the dynamics and functional interactions between regions using a simultaneous recording.

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

We thank Florent Meyniel and Tianming Yang for their comments on the manuscript, and Xinjian Jiang, Jian Jiang, and Juntao Feng for experimental assistance. This work was supported by the National Science and Technology Innovation 2030 Major Program 2021ZD0204204 to WF, the Shanghai Municipal Science and Technology Major Project 2021SHZDZX to LW, the Lingang Laboratory Grant LG202105-02-01 to LW, the Strategic Priority Research Programs XDB32070201, the Strategic Priority Research Programs XDB32070201 to LW, and the National Natural Science Foundation of China 32100830 to WF.

All animal procedures were approved by the Animal Care Committee of Center for Excellence in Brain Science and Intelligence Technology, Institute of Neuroscience, Chinese Academy of Sciences (Permit Number: CEBSIT-2020034).

1. Received: December 6, 2021
2. Preprint posted: December 7, 2021 (view preprint)
3. Accepted: October 23, 2022
4. Accepted Manuscript published: October 24, 2022 (version 1)
5. Version of Record published: November 8, 2022 (version 2)

© 2022, Qi, Fang et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.