Tag: Decision

A Drosophila computational brain model reveals sensorimotor processing

[ad_1]

The Drosophila melanogaster central brain connectome, comprising more than 125,000 neurons and 50 million synaptic connections, allows brainwide analyses of how the fly processes sensory information^{1,2,11,12,13,14,15,16}. Although the Drosophila brain is comparatively small, a single Drosophila neuron may be connected synaptically to hundreds of downstream neurons, and interpreting how this connectivity relates to behaviour is not straightforward.

To model the neural circuit mechanisms that generate behaviour, we implement a simple leaky integrate-and-fire model using the connection weights derived from the entire adult Drosophila central brain connectome of reconstructed electron microscopy neurons^{1,17,18,19,20}, as well as neurotransmitter predictions for each neuron^3,21. In this model, spiking of a neuron alters the membrane potential of downstream neurons in proportion to the connectivity from the upstream neuron¹⁷ (Fig. 1a and Methods); if a downstream neuron’s membrane potential reaches the firing threshold, that neuron, in turn, fires.

Fig. 1: The computational model accurately predicts neurons that respond to sugar stimulation and neurons required for proboscis extension to sugar.

a, Computational model schematic. Activation of the grey neuron at the times indicated by the arrows causes depolarization of the green and purple neurons in proportion to their connectivity from the grey neuron. Upon reaching the firing threshold, a neuron spikes, and its membrane potential is reset. b, Schematic of the proboscis extension response: unilateral sugar presentation causes proboscis extension towards the side of the fly with sugar. Arrow widths are proportional to the number of synapses connecting each set of neurons. c, Predicted MN9 firing rate of either the ipsilateral or contralateral MN9 in response to unilateral left hemisphere sugar GRN activation. Error bars represent s.d., Mann–Whitney U test. d, Heatmap depicting the predicted firing rates in response to unilateral 10–200 Hz sugar GRN firing. The y axis is ordered by firing rate at 200 Hz sugar activation, and depicts the top 200 most active neurons. e, Heatmap depicting the predicted MN9 firing rate when the top 200 responsive neurons are activated at 25–200 Hz; in d and e, the grey squares represent neurons that respond to sugar. f, Heatmap depicting the change in the contralateral MN9 firing rate in response to activation of sugar GRNs, while individually silencing each of the top 200 responsive neurons. The MN9 firing rate for each silenced neuron is normalized relative to the firing of MN9 when no neurons are silenced for each GRN activation frequency; In e and f, the y axis is ordered as in d. g, Histogram of the non-GRNs in f at 50 Hz. h, Venn diagram depicting the intersection between neurons predicted to activate MN9 and neurons predicted to cause a 20% decrease in MN9 firing when silenced.

We implemented this model in the spiking neural network simulator Brian2 (ref. ²²) using all Flywire neurons. The baseline firing of each neuron in our model is 0 Hz. All biophysical parameters are taken from previous D. melanogaster modelling or electrophysiology efforts^18,19, from the synaptic weights from the Flywire connectome^1,21 and from the neurotransmitter predictions³, except for W_syn—the single free parameter of the model, which corresponds to the magnitude of the change in downstream membrane potential that results from a single excitatory or inhibitory synapse (Methods). By driving activity in a sparse set of neurons, the model predicts changes in downstream firing.

We examined the model’s ability to predict circuit activity in two systems: feeding initiation and antennal grooming. We began by examining the Drosophila feeding initiation circuit because it has well-defined taste sensory inputs and motor outputs that are contained in the Flywire electron microscopy volume. Thus, computational modelling of the feeding initiation circuit permits analysis of an entire sensorimotor circuit, in contrast to other sensorimotor circuits that require descending neurons, which are incomplete in the Flywire volume. Furthermore, gustatory sensory neurons that respond to sugar, water and bitter tastes have been identified in the electron microscopy volume²³, permitting a detailed analysis of how these modalities interact. In addition, extensive experimental analysis provides a ground truth for computational studies^{4,5,23,24,25,26,27,28,29,30,31,32}. We further assessed the performance of the model in another well-defined but non-overlapping circuit—the antennal grooming circuit—as an independent evaluation of the model^6,7,8,9,10. As with the feeding initiation circuit, the antennal grooming circuit has well-defined sensory inputs, and a discrete, easily quantified behavioural output: antennal grooming behaviour. In both circuits, we tested specific predictions that the computational model generated using cell-type-specific genetic tools, optogenetics and functional imaging. We find that the model makes predictions consistent with our empirical observations, such as identification of neurons required for behavioural output. Thus, our computational model reduces the vast complexity of the connectome into simple, intuitive circuits.

In Drosophila feeding initiation, detection of appetitive substances in hungry flies results in proboscis extension and consumption³³. Gustatory receptor neurons (GRNs) on the body surface of the fly, including the labellum (tip of the proboscis) or the legs, respond directly to tastants and project to the primary taste centre of the insect brain—the suboesophageal zone (SEZ)^{23,33,34,35,36,37}. GRNs respond to specific taste categories, such as appetitive sugar or aversive bitter compounds, resulting in acceptance (that is, proboscis extension and feeding) or avoidance, respectively^24,33,35,38.To examine the neural circuits that influence feeding in response to taste detection, we focussed on four GRN categories: sugar, water, bitter and a fourth GRN category labelled by the ionotropic receptor Ir94e. Ir94e neurons respond to salt and the presentation of male genitals, but the exact tastants Ir94e neurons respond to are not well understood^31,33,37,38. These GRNs have been identified and classified previously in the electron microscopy brain volume²³; we verify and expand on this classification by clustering on the basis of connectivity and comparing this clustering with response properties of second-order neurons (Extended Data Fig. 1a–c and Methods).

When a fly encounters sugar, activation of appetitive GRNs results in activation of proboscis motor neurons (MNs) (Fig. 1b). The proboscis consists of three segments: the rostrum, the haustellum and the lip-like labella, controlled by the activity of 16 MNs³². We find that computational activation of labellar sugar-sensing GRNs activates several proboscis MNs involved in feeding, including MNs 6, 8, 9 and 11 (Fig. 1c, Extended Data Fig. 1d and Supplementary Table 1)³⁹. Consistent with the model’s predictions, MN9 and MN11 have been shown previously to respond to sugar stimulation in vivo^5,39. In total, we find that the computational model can model a complete sensorimotor transformation.

To assess the ability of our computational model to predict the composition and function of the feeding initiation circuit, we focussed specifically on the activity of MN9, which controls rostrum lifting during proboscis extension^5,32. The rostrum is the largest portion of the proboscis, permitting quantification of MN9 activity by measuring rostrum lifting. Although the exact correlation between MN9 firing rate and rostrum lifting is not known, we assume that increased MN9 firing rates correspond with increased rostrum lifting probability. Remarkably, unilateral sugar GRN activation activates the contralateral MN9 more strongly compared with the ipsilateral MN9 when either the left (Fig. 1c) or the right (Extended Data Fig. 1d) hemisphere GRNs are activated, consistent with behavioural experiments showing that unilateral taste detection on the legs promotes proboscis extension that is curved and directed towards the food source^26,40. Thus, we show that in silico sensory activation produces MN activity that is consistent with the observed behaviour of the fly taste sensorimotor circuit.

To confirm that our computational activation of MN9 depends on the actual connectivity weights determined from the fly connectome, we tested whether distorting synaptic weights would impair the ability of sugar sensory neurons to activate MN9. In these experiments, connectivity weights were shuffled randomly (while maintaining the global connectivity weight distribution). Although modelling using the correct connectome results in robust activation of MN9 in 100% of simulations when sugar-sensing neurons are activated at 100 Hz, only 1 of 100 shuffled simulations did (Supplementary Table 1d). Therefore, the predictive accuracy of our computational model depends on the actual connectivity weights of the fly connectome.

We next examined whether the computational model could accurately predict the neuronal cell types that are known to compose the feeding initiation circuit⁴. We first examined the neural network activated upon unilateral sugar GRN activation. We note that, given the variety of assumptions the model relies upon, absolute firing rate predictions are unlikely to be accurate; therefore, we examined network activity upon sugar GRN activation ranging from 10 to 200 Hz (Fig. 1d). We find that increasing sugar GRN firing rate increases activity of MN9, as well as MNs 6, 8 and 11. Of the 127,400 neurons modelled, we found that 45 are predicted to respond to 10 Hz sugar GRN activation, and 455 to 200 Hz (Supplementary Table 1). Activated neurons are defined as neurons that have greater than 0 Hz firing. Thus, the computational model predicts a large network activated by sugar taste detection that includes known sugar-responsive MNs.

Sugar taste detection influences activity in nutritive state and memory circuits, and modulates a broad range of behaviours, including feeding, oviposition and foraging^33,37,41. To specifically evaluate the subset of predicted sugar-responsive neurons that influence feeding initiation, we performed two further in silico experiments. First, as a strategy to identify neurons that drive feeding initiation, we computationally stimulated each of the top sugar-responding neurons in the network to identify those that drive activity in MN9 (Fig. 1e). Second, to identify neurons required for feeding initiation to sugars, we computationally activated sugar GRNs, silenced each of the top 200 sugar-responsive neurons one at a time, and measured the change in predicted MN9 firing (Fig. 1f,g). For these silencing experiments, we activated sugar-sensing neurons at frequencies ranging from 50 to 120 Hz in 10 Hz increments. Neurons that our model predicts to be required for feeding initiation will have decreased MN9 firing when silenced. We defined neurons predicted to cause a silencing phenotype as any neuron whose silencing causes MN9 firing to be 80% or lower compared with control MN9 firing at any of the eight sugar activation frequencies tested (50, 60, 70, …120 Hz; Fig. 1f). In general, silencing of individual neurons had the greatest effect when sugar GRNs were activated at low frequencies, implying greater redundancy in the circuit as sensory stimulation increases. In total, our analyses identified 47 neurons predicted to be sugar-responsive, and sufficient for feeding initiation. Of these 47 neurons, 14 are also predicted to be required for MN9 activity (Fig. 1h).

We next evaluated whether the predicted neurons for feeding initiation include neurons shown experimentally to participate in feeding initiation behaviour. Previous experimental studies identified ten neural classes that respond to sugar, and are sufficient for proboscis extension⁴ (Fig. 1b and Extended Data Fig. 2a–c). Our computational model correctly predicts that all ten cell types respond to sugar (Supplementary Table 2). Of these ten neurons, eight are predicted correctly to be sufficient to activate MN9 (ref. ⁴) (Supplementary Table 2). We previously found that five of the ten are required for sugar feeding initiation⁴ (Supplementary Table 2). Of these five, three are predicted by our computational model to cause a greater than 20% decrease in MN9 firing, and one of the others is predicted to cause a statistically significant decrease in MN9 firing, but less than 20%, when silenced. Although the model predictions generally match previous experimental results, there are some deviations. For example, the model fails to correctly predict that the Phantom cell type will activate MN9 (ref. ⁴) (Fig. 1e, Extended Data Fig. 2c and Supplementary Table 2). This cell type is predicted to be inhibitory. Phantom strongly synapses onto Scapula—a neuron that is also predicted to be inhibitory; Scapula, in turn, synapses onto Roundup, the pre-MN with the strongest predicted silencing phenotype. We speculate that activation of Phantom inhibits Scapula, potentially permitting Roundup and MN9 firing. Because the basal firing rate of all neurons in the model is 0, activation of inhibitory neurons in the model, in the absence of other input, cannot alter the firing of downstream neurons. A further explanation for incorrect predictions could be neuromodulation, which is not accounted for in our model. Particular neurons may be subject to neuromodulation, causing their activity to be different from predictions on the basis of connectivity. Alternatively, neurons that express neuromodulators may be poorly modelled. We speculated that the Usnea cell type, which has a strong experimental activation and silencing phenotype⁴ yet is not predicted to be either necessary or sufficient for proboscis extension, might be neuropeptidergic. To test this, we performed cell-specific knockdown of the gene Amontillado—a prohormone convertase required for neuropeptide processing in Drosophila^42,43. Knockdown of Amontillado phenocopied the Usnea silencing phenotype (Extended Data Fig. 1e), indicating that Usnea activity may require neuropeptide processing. Additionally, incorrect neurotransmitter predictions or other assumptions of the model may explain discrepancies between the prediction of our model and our experimental results. Despite these limitations, overall, this analysis demonstrates that our computational model correctly identifies known neurons in a sensorimotor circuit.

As an independent assessment of whether the computational model accurately predicts neurons that elicit MN9 activity, the output of our sensorimotor circuit, we compared optogenetic activation phenotypes with their corresponding computational activation phenotypes. To do this in a non-biased way, we performed a screen in which we optogenetically activated individual neuronal cell types with split-GAL4 lines and monitored the activity of MN9. The SEZ split-GAL4 collection labels 138 cell types in the SEZ—the primary feeding region of the brain⁴⁴. We identified 106 of these labelled cell types in the Flywire volume. Next, we crossed these split-GAL4 lines to create flies expressing the light-gated cation channel CsChrimson. We then optogenetically activated these neurons, and measured whether MN9 is activated by observing rostrum extension. We compared the predicted in silico MN9 activation phenotypes of these cell types when we activate them between 10 Hz and 200 Hz with the actual optogenetic activation MN9 phenotypes we observed. When we activate each cell type at 50 Hz, 11 are predicted to activate MN9 (that is, elicit MN9 firing greater than 0 Hz); notably, 10 of 11 of these cell types actually do elicit rostrum extension when optogenetically activated (Fig. 2a–c and Supplementary Table 3). Furthermore, of the 95 predicted not to elicit proboscis extension due to 50 Hz activation, just 4 have non-zero rostrum extension. Activation of these cell types at 200 Hz, rather than 50 Hz, results in the addition of five false positives. At 10 Hz activation, six cell types are predicted to cause MN9 activation; of these five, six do indeed cause proboscis extension. Thus, the computational model can predict the activation phenotypes of a non-biased sample of cell types at greater than 90% accuracy.

Fig. 2: The computational model predicts neurons that cause proboscis extension.

a, Predicted MN9 firing rates when each of 106 cell types are activated computationally at 50 Hz. Cell types are ordered by predicted MN9 firing rate. b, Fraction of flies extending the rostrum—the segment of the proboscis controlled by MN9—in response to optogenetic activation; cell types are ordered as in a. Where several split-GAL4 lines for a cell type were tested, the line with the highest extension rate is plotted; n = 10 flies per cell type. c, Confusion matrix showing the accuracy of MN9 activation predictions and the number of cell types in each category. The rostrum was predicted to extend if MN9 was predicted to have non-zero firing as a result of 50 Hz cell-type activation.

The accuracy of the model indicates that it provides a powerful platform to discover how different taste modalities are processed to influence feeding initiation. We first tested whether the model can predict the response to coactivation of both an attractive sugar stimulus and an aversive bitter stimulus. Bitter detection inhibits proboscis extension motor activity⁴⁵ (Fig. 3a). Indeed, the addition of bitter GRN activity to sugar GRN activation in our model resulted in an inhibition of MN6 and MN9 (Fig. 3b and Supplementary Table 4). We previously found, using calcium imaging, that bitter GRN activation inhibits the sugar pathway at the level of pre-MNs⁴, consistent with the predictions of the computational model (Supplementary Table 4).

Fig. 3: The computational model correctly predicts that Ir94e neurons are aversive but fail to inhibit proboscis extension to a strong sugar stimulus.

a, Schematic outlining the previously known and unknown roles of sugar, bitter and Ir94e neurons. Question marks indicate that the exact substrate(s) that activate Ir94e neurons are not known, nor is it known whether Ir94e activation influences proboscis extension. b,c, Heatmap depicting the predicted MN9 firing rates in response to the combination of sugar GRN firing and bitter (b) or Ir94e (c) GRN activation. d,e, Fraction of flies exhibiting PER upon 50 mM sucrose stimulation or 1 M sucrose stimulation when Gr66a/bitter GRNs (d) or Ir94e GRNs (e) are optogenetically activated. Red bars indicate red light condition; n = 26–32, see Supplementary Table 9 for exact values. Mean ± 95% confidence intervals using Wilson’s score interval, Fisher’s exact test. f, Venn diagram showing the number and overlap of neurons that respond to sugar GRN (green) or water GRN (blue) activation that elicits 40 Hz MN9 firing, as well as bitter GRN (red) or Ir94e GRN activation (purple) activated to reduce 40 Hz MN9 firing to 1 Hz.

We next examined the predicted circuit activity caused by GRNs labelled by the ionotropic receptor Ir94e; these neurons have been identified previously in the electron microscopy volume²³. Ir94e neurons respond to low salt concentrations and the presentation of male genitals, among other substances^31,38, and are suggested to play a role in mediating attraction to low salt³¹. However, the role they play in proboscis extension has not been described. Notably, the computational model predicted that activation of Ir94e GRNs, rather than promoting MN9 firing, inhibits MN9 firing (Fig. 3c). Therefore, we tested whether optogenetic activation of Ir94e GRNs is sufficient to inhibit proboscis extension, similar to bitter activation. Indeed, we found that optogenetic activation of Ir94e GRNs or bitter GRNs was sufficient to inhibit the proboscis extension to 50 mM sucrose, as our modelling predicted (Fig. 3d,e). Interestingly, we noted a quantitative difference between the model’s predictions for bitter versus Ir94e activation. Strong bitter activation is predicted to eliminate MN9 firing to strong sugar stimulation, but strong activation of Ir94e neurons is not predicted to do so (Fig. 3b,c). We therefore tested the proboscis extension response (PER) to 1 M sucrose while optogenetically activating bitter or Ir94e GRNs. Optogenetic bitter activation eliminated consumption of 1 M sucrose (Fig. 3d), but Ir94e activation did not (Fig. 3e). Thus, we conclude that Ir94e GRN activity inhibits proboscis extension, but fails to fully inhibit proboscis extension to strong sugar stimuli. These results indicate that our computational model can predict previously unknown circuit functions and properties.

Finally, we sought to predict how water taste detection influences feeding initiation. The degree to which sugar GRNs and water GRNs activate pathways that are distinct or shared is unknown. We found that activation of water GRNs in our model activates many downstream neurons that are also activated by sugar stimulation. (Fig. 3f, Fig. 4a and Supplementary Tables 1 and 6). In particular, comparing neurons activated by sugar GRNs with those activated by water GRNs, at a stimulation frequency at which each pathway activates MN9 at 40 Hz, predicted that the sugar pathway activates 377 neurons, while the water pathway activates 391 neurons. Of these, more than half (250), are shared between the two circuits (Fig. 3f and Supplementary Table 4). We also examined bitter responsive neurons and Ir94e responsive neurons at the minimal activation sufficient to reduce 40 Hz MN9 firing to 1 Hz. Only two neurons were common between sugar and bitter activation, and 30 between sugar and Ir94e activation, demonstrating segregation of neurons activated by aversive and appetitive taste (Fig. 3f). This prediction is consistent with our previous calcium imaging experiments demonstrating that, across nine sugar-responsive cell types, zero respond to a mixture of bitter compounds⁴. In contrast, our model predicts central neurons that respond to both sugar and water taste activation, as well as sugar-specific and water-specific neurons, consistent with brainwide calcium imaging studies^29,46.

Fig. 4: The computational model correctly predicts that the sugar and water pathways share components and additively promote proboscis extension.

a, Heatmap depicting the predicting firing rates in response to 20 to 260 Hz water GRN firing. The y axis is ordered by firing rate at 260 Hz water activation. b, Heatmap depicting the predicted MN9 firing rate when the top 200 responsive neurons are activated at 25–200 Hz. c, Heatmap depicting the change in MN9 firing rate in response to activation of water GRNs at the specified firing rate, while individually silencing each of the top 200 responsive neurons. d, The fraction of flies exhibiting PER upon water stimulation. Green bars indicate green light condition; n = 30–50; see Supplementary Table 9 for exact values. Open and filled circles represent whether the computational model predicted a greater than 20% decrease in MN9 firing at 160 Hz water GRN stimulation. e, Heatmap depicting the predicted MN9 firing rates in response to the combination of sugar and water GRN activity. f, The fraction of flies exhibiting PER upon water stimulation. n = 39–40. d,f, Mean ± 95% confidence intervals, Fisher’s exact test.

To identify interneurons that compose the water feeding initiation circuit, we used the computational model to analyze the water-responsive neurons that influence MN9 activity. We stimulated the top 200 neurons that are predicted to respond to water, and identified the subset that computationally activates MN9 (Fig. 4b). Next, we also computationally activated water-sensing GRNs, silenced each water-responsive neuron, and monitored the change in MN9 activity (Fig. 4c and Extended Data Fig. 2d). Using our computational model, we identified 39 water-responsive neurons that are also sufficient for MN9 activation. Of these 39, 30 are also predicted to be activated by sugar GRNs (Extended Data Fig. 2d). Furthermore, we identify nine neurons predicted to be both necessary and sufficient for water feeding initiation (Extended Data Fig. 2e). As with sugar, we defined a neuron predicted to be required for water feeding initiation as any neuron that, when silenced, caused MN9 firing to be less than 80% of that of the unsilenced control.

To test these predictions experimentally, we performed calcium imaging on two neurons predicted to respond to water: Fudog and Zorro. We found that both neurons indeed responded to water (Extended Data Figs. 1c and 3a). Additionally, we examined six neurons predicted to have water silencing phenotypes. Five of these, when silenced optogenetically, indeed decreased significantly proboscis extension to water, while a sixth, G2N-1, did not (Fig. 4d and Extended Data Fig. 3b,c). We also examined five neurons that respond to computational water activation, but are not predicted to cause a water silencing phenotype. Of these five neurons, four did not have a water silencing phenotype, as predicted, although one, Usnea, did decrease proboscis extension significantly when silenced with GtACR1 (Fig. 4d and Extended Data Fig. 3b,c).

Our computational model predicts that the water and sugar pathways share a common set of neurons (Extended Data Fig. 2d,e). Do these shared neurons contribute to feeding initiation? Our calcium imaging experiments (Extended Data Figs. 1c and 3a) combined with previous experiments⁴ confirm that five neurons predicted to respond to sugar and water do respond to both sugar and water in vivo: Clavicle, Fudog, Phantom, Rattle and Zorro. Moreover, four of these neurons had been shown previously to be sufficient for proboscis extension, and three are also required for sugar feeding initiation⁴. All three are among the neurons we found experimentally to be required for feeding initiation to water, as predicted (Fig. 4d). Furthermore, the two other cell types we found experimentally to be required for water—Bract and Roundup—are also predicted to respond to both water and sugar (Supplementary Tables 1 and 6), and have been found to respond to sugar⁴. However, previous calcium imaging studies did not identify water responses in these two cell types⁴. This discrepancy may reflect the greater sensitivity of the behavioural silencing experiments compared with calcium imaging of water responses⁴. Finally, a further cell type, Usnea, has been shown to respond to water, but not sugar⁴; our model correctly predicts Usnea responds to water, but incorrectly predicts that it will also respond to sugar. Usnea has previously been shown to be required for feeding initiation to sugar, and our Amontillado RNAi experiments indicate that it may be neuropeptidergic. We find that it is also required for proboscis extension to water (Fig. 4d). Usnea synapses directly onto both sugar and water GRNs (Extended Data Fig. 1b), and may tune the response of these neurons. Thus, we identify a set of neurons involved in the processing of both sugar and water.

To explore the relationship between the water and sugar pathways, we computationally activated both sugar and water GRNs simultaneously and examined the effect on MN9. Our computational modelling predicts that activation of water and sugar GRNs work synergistically to promote MN9 firing (Fig. 4e and Extended Data Fig. 3d). If sugar and water do act synergistically, then both sugar GRNs and water GRNs may be involved in water consumption. Only water GRNs have been implicated in proboscis extension to water; we asked whether sugar GRNs might also be required. Indeed, silencing of sugar GRNs reduced the fraction of flies that extended their proboscis to water (Fig. 4f and Extended Data Fig. 3e). Further, silencing water-sensing neurons reduced consumption of 50 mM sucrose, although a confound is that these water-sensing neurons are known to respond to this concentration of sucrose (Extended Data Fig. 3f). In total, our computational modelling, optogenetic behaviour experiments and functional imaging indicate that the water and sugar pathways share, at least in part, common components to form an appetitive consumption pathway.

To test the general applicability of the computational model to study sensorimotor processing, we sought to determine whether it could predict circuit properties in another system—the well-studied antennal grooming circuit^6,7,8,9,10. In this system, activation of a set of mechanosensory neurons in the Johnston’s organ—a chordotonal organ in the antennae—elicits grooming of the antennae^8,47 (Fig. 5a). These mechanosensory neurons, abbreviated JONs, synapse onto two interneuron types, named antennal grooming brain interneurons 1 and 2 (aBN1 and aBN2), which in turn synapse onto two descending neurons, aDN1 and aDN2 (ref. ⁸). There is a single aBN1 per hemisphere, while there are several aBN2 neurons per hemisphere. Each of these cell types—aBN1, aBN2, aDN1 and aDN2—are sufficient for antennal grooming, while aBN1 and aBN2 are each at least partially required for antennal grooming⁸.

Fig. 5: The computational model correctly identifies key neurons in the antennal grooming circuit as well as subtype circuit responses.

a, Schematic of the antennal grooming circuit. Arrows represent known functional connectivity⁸. Grey oval around aDNs indicates that JONs activate aDNs, but exactly which aDNs are not known. b, Heatmap depicting the predicting firing rates in response to 20–220 Hz JON firing; 147 JONs were activated, and are the neurons that have the highest firing rates. Neurons are ordered by firing rate at 220 Hz. c, Heatmap depicting the predicted aDN1 firing rate when the top 300 responsive neurons are activated at 25–200 Hz. d, Heatmap depicting the change in aDN1 firing rate in response to activation of JOs at the specified firing rate, while individually silencing each of the top responsive neurons. e, Histogram of the predicted change in aDN1 firing rate as a result of silencing each non-JONs, when JONs are activated at 140 Hz. The y axis depicts the number of neurons in each bin. Neurons previously identified are labelled. f, Venn diagram depicting the overlap between neurons predicted to be sufficient to activate aDN1 at greater than 2 Hz and neurons required for aDN1 activation. g, JO subtype connectivity onto aBN1 and predicted aBN1 firing in response to JO activation at the specified rate. Error bars, s.d. h, Calcium imaging of aBN1 in response to optogenetic activation of each subtype. The ΔF/F average ± s.e.m. is shown; n ≥ 5 flies tested. Shaded bar indicates when a red light pulse was delivered.

We first sought to test whether the computational model could identify the previously described neurons in the circuit. We activated a set of 147 previously identified JONs of the JO-C, JO-E, JO-F and JO-m subclasses^8,47. Indeed, the model identified that aBN1, aBN2, aDN1 and aDN2 respond to JON activation (Fig. 5b and Supplementary Table 5). To determine which of these JON-responsive neurons might drive antennal grooming, we computationally activated these neurons and asked whether they could elicit activity in either of the two descending neurons that evoke antennal grooming: aDN1 or aDN2 (Fig. 5c and Extended Data Fig. 4a). Next, we asked, among the top neurons predicted to respond to JON activation, which are required for activation of aDN1 or aDN2 (Fig. 5d and Extended Data Fig. 4b). Notably, only four neurons, beyond aDN1 itself, were identified that could elicit aDN1 activity: aBN1, aDN2 and two other neurons that elicited less than 2 Hz aDN1 activity (Fig. 5c and Supplementary Table 7). Moreover, only three neurons, besides aDN1 itself, were identified that reduced aDN1 activity by more than 20% at 140 Hz JON activation: aBN1; a descending member of the BN2 class; and aDN2 (Fig. 5d–f). Thus, the computational model identifies members of each of the previously identified critical nodes of the antennal grooming circuit purely from knowledge of the sensory inputs and descending outputs.

We next tested how different JON subpopulations influence antennal grooming. JONs send their projections to the antennal mechanosensory and motor center in the ventral brain. JO-C and JO-E neurons respond to antennal vibrations and project medially into the antennal mechanosensory and motor center, while JO-F neurons project into a distinct region⁸. Optogenetic activation of both JO-CE and JO-F neurons is sufficient to trigger antennal grooming, but it is not known whether these two populations generate distinct patterns of downstream firing. Both JO-CE and JO-F neurons synapse onto aBN1 (103 and 78 synapses, respectively; Fig. 5g), raising the possibility that they elicit grooming by activating aBN1.

Our computational model predicts that, whereas JO-CE neurons will elicit robust aBN1 activity, JO-F neurons will not, despite synapsing directly onto aBN1 (Fig. 5g). To test this prediction, we optogenetically activated each population of JONs and performed calcium imaging in aBN1. Consistent with the prediction of this model, JO-CE activated aBN1 robustly, but JO-F neurons did not (Fig. 5h). Why do JO-F neurons fail to activate aBN1 robustly? We identified three putative inhibitory neurons that are directly postsynaptic to JO-F neurons and synapse directly onto aBN1. Computational silencing of these three neurons permits JO-F neurons to activate aBN1, but this remains to be tested empirically (Extended Data Fig. 4c,d). Our analysis of the antennal grooming circuit demonstrates that our computational model can provide insights into complex circuits, purely from knowledge of sensory input and descending output. We demonstrate that modelling brain circuits purely from connectivity and neurotransmitter identity is sufficient to reliably describe, at least at a coarse level, entire sensorimotor transformations.

In conclusion, we report a computational model on the basis of connectivity and neurotransmitter predictions of the entire fly connectome that can predict circuit neural activity, the neurons required for activation of output neurons and the integration of several sensory modalities. We use the model to create predictions of the sugar, water, bitter and Ir94e pathways and validate many of these predictions experimentally. We show that the Ir94e neurons, previously considered to be attractive, instead inhibit proboscis extension. The results of our modelling indicate that sugar, bitter and Ir94e GRNs activate generally distinct populations of neurons. In contrast, sugar and water GRNs activate many of the same central neurons as well as sugar-specific and water-specific neurons. In addition, we recapitulate the antennal grooming circuit purely from sensory input and descending output, and identify a subpopulation of JONs that, despite strong connectivity onto aBN1, fail to activate it. These studies demonstrate the power of computational modelling to explain sensory processing features in complex networks.

Our analysis of the taste and antennal grooming circuits shows we can model local sensorimotor transformations in the taste and antennal grooming circuits. The computational model, implemented in the widely used Brian2 library²², allows for perturbations that are easily interpreted. We believe our computational model will be a useful tool for the study of sensorimotor transformations and the exploration of interactions between overlapping neural pathways (for example, sweet-bitter, sweet-water and so on).

[ad_2]

Source link

October 2, 2024
Brain-wide dynamics linking sensation to action during decision-making
[ad_1]
Animals

All experiments were performed under the UK Animals (Scientific Procedures) Act of 1986 (PPL: PD867676F) following local ethical approval by the Sainsbury Wellcome Centre Animal Welfare Ethical Review Body. A total of 21 C57BL/6 J male mice (age = 34.5 ± 15.8 weeks (mean ± s.d.)) were used for electrophysiological recordings. Fifteen mice first underwent head-fixed behavioural training prior to acute electrophysiological recordings (see ‘Task and training stages’), and six mice (untrained mice) only underwent habituation to the recording setup prior to acute electrophysiological recording.

Prior to behavioural training and recordings, all mice were implanted with a head-fixation bar under approximately 1.5% isoflurane and administration of Meloxicam (5 mg kg⁻¹) to allow for head-fixation during behavioural training and electrophysiological recordings.

During training, mice were co-housed with littermates in individually vented cages. After implantation of the recording chamber, mice were singly housed to protect the implant. Mice were housed in reversed day–night cycle lighting conditions, with the ambient temperature and humidity set to 23 °C and 56% relative humidity, respectively.

Behavioural task

The design of the behavioural task was as previously described in ref. ¹⁴. In brief, mice were head-fixed and placed on a polystyrene wheel. Two monitors (21.5 inch, 1,920 × 1,080, 60 Hz) were placed on each side of the mouse at approximately 20 cm from the mouse head. The monitors were gamma corrected to 40 cd m⁻² of maximum luminance using custom MATLAB scripts utilizing PsychToolbox-3. The stimulus presentation was controlled by custom written software in MATLAB utilizing PsychToolbox-3. The visual stimulus was a sinusoidal grating with the spatial frequency of 0.04 cycles per degree resulting in 3 grating periods shown on a screen. Each trial began with a presentation of a grey texture covering both screens. After a randomized delay (at least 3 s plus a random sample from an exponential distribution with the mean of 0.5 s), the baseline stimulus appeared. The TF of the grating was drawn every 50 ms (3 monitor frames) from a lognormal distribution, such that log₂-transformed TF had the mean of 0 and s.d. of 0.25 octaves and the geometric mean of 1 Hz. The direction of drift was randomized trial to trial between upward or downward drift. The sustained increase in TF, referred to in the text as change period, occurred after a randomized delay (3–15.5 s) from the start of baseline period and lasted for 2.15 s. For early and late blocks training (stage 8), change period times were sampled between [3, 8] s and [10.5, 15.5] s, respectively, with the delay from the earliest allowed change period sampled from an exponential distribution with a mean of 4 s. Random 15% of trials were assigned as no-change trials and did not have a change period. For stage 8 training, 10% of trials were designated to be probe trials and had a change time drown from the distribution of the other block type. Because there were no qualitative differences in neural TF pulse response between early and late blocks (data not shown) we have combined data from both block types for analyses throughout this manuscript. Findings related to stage 8 (early and late blocks) will be presented in an upcoming paper.

Mice were trained to report sustained increases in TF by licking the spout to trigger reward delivery (drop of soy milk). Licks that occurred outside of the change period are referred in the text as early licks. If mice moved on the wheel (movement exceeding 2.5 mm in a 50-ms window) in either direction, the trial was aborted (stages 7 and 8). If mice did not lick within 2.15 s from the change onset, the trial was considered a miss trial.

Training stages

Following the implantation of the headplate, mice were allowed to recover for a week. After that, mice went through several stages of training:
1. (1)
  
  Mice were handled for 3 to 7 days, until mice were comfortable with being handled by the experimenter. During this stage mice were also habituated to being restrained by being placed into a soft cloth for a short period of time. After the brief restraints they were given a small amount of soy milk as reward.
2. (2)
  
  Next, mice were put on food restriction. Mouse weight was monitored daily with the amount of food given adjusted per mouse to keep them sufficiently motivated for getting rewards and keep their weight no lower than 85% of the original weight prior to food restriction.
3. (3)
  
  Next, mice were head-fixed and placed on the running wheel of the behavioural training setup with the monitors turned off. Mice were allowed to freely run on the wheel, but not encouraged to. Typically, there were 3 habituation sessions, with the duration progressive increasing from 15 to 45 min.
4. (4)
  
  Next, mice were introduced to the visual stimuli used in the task. Mice were initially shown only trials with two largest changes of TF (2 and 4 Hz, lasting 2.15 s), followed by a reward auto-delivery 1.5 s after the change onset. After mice started to robustly make licks during the change period that preceded the reward auto-delivery, they were transitioned to the next stage.
5. (5)
  
  Here only hit trials were rewarded, early licks and running did not result in termination of the trial.
6. (6)
  
  After mice robustly detected strong changes in the previous step, we introduced trials with weaker changes in TF (1.25 Hz, 1.35 Hz and 1.5 Hz). Additionally, a consequence of an early lick outside of the change period was a mild air-puff to the mouse’s right cheek and a termination of the trial.
7. (7)
  
  After mice detected weaker changes as well (assessed as higher hit rate compared to no-change trials), they were transitioned to the next stage where in order to initiate the trial start (start of the baseline stimulus), mice were required to remain stationary on the running wheel for at least 3 s plus a random sample from an exponential distribution with the mean of 0.5 s. Additionally, after the trial start, a trial was aborted as a consequence of a movement on the wheel.
8. (8)
  
  Finally, after mice reached sufficient proficiency at the previous stage, early and late blocks were introduced. During the session start, a block type was randomly chosen. A block was defined as a period of the session during which a mouse completed 30 hit trials. After completion of a block of trials, the block type was switched to the other block type (early to late or vice versa).
Six mice that were used in the untrained control experiment (Fig. 4e–h) went through training stages 1–3 above. Following that, they were shown the same stimuli as the trained mice, with the difference that their movements on the wheel or licking the spout did not terminate a trial nor trigger reward. Instead, they were given rewards at random times with inter-reward intervals drawn from the uniform distribution of 60 ± 15 s.

Behavioural setup and data acquisition

Reward delivery (soya milk) was controlled by a solenoid pinch valve (161P011, NResearch) and delivered to the mouse via a spout positioned in front of it. Mouse licking the spout was measured by a piezo element (TDK PS1550L40N) coupled to the spout and amplified with a custom-made amplifier system. Running wheel movement was measured with a rotary encoder (model Kübler) that was connected to the wheel axle. All behavioural data and events, such as piezo signal voltage trace, valve or change period on/off state, etc., were acquired via analogue and digital channels of PXIe-6341 acquisition card (National Instruments) with SpikeGLX (https://github.com/billkarsh/SpikeGLX) at 8,474 Hz.

Behavioural data analysis

Psychometric performance, reaction times and lick-triggered stimulus average

Psychometric curves were calculated per session by counting the amount of hits relative to all trials where mice did not early lick nor abort. Mean hit rates (performance) and parametric 95% confidence intervals (s.e.m. × 1.96) of hit rates were calculated across sessions (n = 114) per change size. Mean reaction times and parametric 95% confidence intervals were calculated across sessions (n = 114) per change size, and p-values were estimated from t-tests.

Lick-triggered stimulus average was estimated by extracting the TF pulses from −1.5 to 0 s preceding early licks and averaged across all trials, revealing mean stimulus TF prior to early licks. Parametric 95% confidence intervals were estimated by calculating the s.e.m. of TF values at each 50 ms bin (TF pulse resolution) prior to an early lick and multiplying the s.e.m. by 1.96.

Simple behavioural leaky-integrator model

In order to formally test if mice behaviourally integrated stimulus evidence (TF pulses) over time in our task, we constructed a simple behavioural leaky-integrator model with two adjustable parameters: decay time (τ), and threshold. We fitted these two parameters by estimating which decay time and threshold predicted most early lick times (from 2 s after trial start, to exclude trial onset licks) correctly for each mouse and then determined the average best-fit decay time and threshold values across mice. For each early lick trial, we calculated the integrated log-scaled TF with decay across the entire trial up until the early lick.

For each early lick trial, we then estimated whether a threshold crossing of the integrated TF had been predicted within a second preceding an actual early lick onset. If this was the case, we considered the model to have predicted the early lick time. If not, we considered the model to not have correctly predicted that trial. We did this for all early lick trials, using a 58 × 151 parameter space: 58 possible decay times spanning from 0.05 s decay time (that is, no integration) to 1,000 s decay time (that is, perfect integration): (50 log-spaced decay times spanning 0.050–3 s, as well as 8 additional very long decay times: 4, 5, 6, 7, 8, 9, 20, 1,000 s), and 151 linearly spaced thresholds spanning [0.01–0.16]. Significance testing of best decay time across mice (that is, larger than no integration (0.05 s)) was done with a t-test.

We also tested if the best-fit decay/integration time parameter estimated from predicting early lick times also outperformed a model with no integration when predicting single-trial hit reaction times (that is, a trial type which the parameters were not optimized on). We did this by comparing actual and predicted reaction times per change size, and calculated Pearson’s correlation between actual reaction times and predicted reaction times per change size. We calculated this by either looking at all reaction times, or only including a subset of trials with reaction times under a defined value (that is, reaction-time cut-off). This was done to better detect if any of the models specifically struggled to predict very late reaction times which may be modulated by non-sensory factors such as such as inattention or lack of engagement. Finally, for significance testing (that is, paired t-test) of whether a model with no integration (decay time = 0.05 s) versus a model with the best-fit decay/integration time (estimated from early lick trials as described above), were significantly different at predicting single-trial reaction times, we z-scored actual and predicted reaction times per change size (to account for change size mean reaction-time differences), and calculated the correlation between all actual reaction times (1 s reaction-time cut-off) and all predicted reaction times of a model with or without integration per mouse, and performed a paired t-test (across mice) of the correlation values from integration versus no integration models.

Outlier detection agent

To test whether mice accumulate evidence over time or merely respond to the instantaneous stimulus, we formulated a null model where behavioural responses are produced via a stochastic outlier detection strategy. Here, an internal decision occurs when a noisy sensory representation of the stimulus crosses a decision boundary, and a response occurs after a stochastic delay. The response is triggered by a single, instantaneous value of the stimulus. However, owing to the stochastic delay, responses may show a gradually decaying statistical dependence on the stimulus history, and may even mimic evidence accumulation strategies such as integration⁴².

Model. According to the outlier detection model, behavioural responses are generated independently for each trial as follows. Let s_i be the stimulus amplitude (log TF) at each time point t_i. We chose time points to correspond with video frames of the stimulus, which were presented at 60 Hz (3 frames per TF pulse). At each time point, a noisy sensory representation Z_i is formed as the sum of the stimulus amplitude and independent and identically distributed (i.i.d.) Gaussian sensory noise ε_i (with mean zero and variance σ²):

$${Z}_{i}={s}_{i}+{\varepsilon }_{i}$$

$${\varepsilon }_{i}\mathop{ \sim }\limits_{\text{i.i.d.}}{\mathscr{N}}\left(0,{\sigma }^{2}\right)$$

An internal decision to respond occurs at time D, given by the first time point where the sensory representation exceeds a decision bound b (or ∞ if the bound is not crossed before the stimulus ends):

$$D=\,\min \{{t}_{i}| {Z}_{i} > b\}\cup \{\infty \}$$

The hazard function of the decision time is thus:

$${H}_{D}\left(d\right)=\prod _{i{\rm{| }}{t}_{i}\le d}p\left({Z}_{i}\le b\right)=\prod _{i{\rm{| }}{t}_{i}\le d}\varPhi \left(\frac{b-{s}_{i}}{\sigma }\right)$$

where Φ is the standard normal cumulative density function (CDF).

A motor response begins at time R, given by the decision time plus an independent, nonnegative stochastic delay ∆ representing the duration of nondecision processes (for example, decision to motor delays):

The delay has a shifted log-logistic distribution with location α, scale β and shape γ, and can be obtained by exponentiating a logistic random variable and then adding a constant. We constrained the location (α > 0) and shape (γ > 1) to give the distribution nonnegative support and a bump-like density that decreases on both sides of the mode. The delay time probability density function (PDF) and CDF are:

$${p}_{\varDelta }\left(\varDelta \right)=\frac{\gamma {\left(\frac{\varDelta -\alpha }{\beta }\right)}^{\gamma -1}}{\beta {\left(1+{\left(\frac{\varDelta -\alpha }{\beta }\right)}^{\gamma }\right)}^{2}}$$

$${F}_{\varDelta }\left(\varDelta \right)=\frac{1}{1+{\left(\frac{\varDelta -\alpha }{\beta }\right)}^{-\gamma }}$$

Because the decision and delay times are independent, the marginal response time distribution is given by the convolution of the decision and delay time distributions. The marginal PDF and CDF of the response time are:

$${p}_{R}\left(r\right)=\sum _{d}{p}_{D}\left(d\right){p}_{\varDelta }\left(r-d\right)$$

$${F}_{R}\left(r\right)=\sum _{d}{p}_{D}\left(d\right){F}_{\varDelta }\left(r-d\right)$$

where the decision time probability mass function (PMF) p_D can be computed from the hazard function H_D above. Because delays are nonnegative, p_Δ(r − d) = F_Δ(r − d) = 0 for all d > r, so the above sums need only be computed over time steps up to the given response time.

The outlier detection model was implemented using custom Python software using the NumPy, SciPy, and PyTorch libraries. All computations involving probabilities were performed in log space, using functions designed to avoid numerical under/overflow.

Fitting. A separate model was fit for each mouse in two stages. We first fit the delay time distribution using only trials with the largest change magnitude, then fit the remaining decision parameters using the entire dataset (excluding the abort trials). This two-stage approach relies on the assumption that delays are identically distributed across trials. In return, it allows more direct estimation of the delay time distribution, providing better ability to distinguish between outlier detection and longer-timescale strategies such as integration.

For each trial i, let n⁽ⁱ⁾ be the number of time points, ${s}^{(i)}=\{{s}_{1}^{(i)},\ldots ,{s}_{{n}^{(i)}}^{(i)}\}$ be the stimulus amplitudes, and c⁽ⁱ⁾ be the time of the change point. For trials where a response occurred, let r⁽ⁱ⁾ be the response time, measured as the onset of facial movement (see ‘Motion onset time estimation’ section) and ℓ⁽ⁱ⁾ be the subsequent lick time (measured at the reward spout).

Fitting the delay time distribution. We assumed that the greatest change magnitude (geometric mean TF 4 Hz) was large enough to trigger an immediate decision at or near the change point. Under this assumption, the delay time on large-change hit trials can be approximated by the reaction time, which can be directly measured as the time elapsed between the change point and the onset of facial movement. Thus, we fit the delay time distribution (shifted log-logistic distribution) to reaction times on large-change hit trials (denoted ${{\mathcal{T}}}_{{\rm{bighit}}}$) by maximum likelihood, subject to the constraints described above:

$$\mathop{\max }\limits_{\alpha > 0,\beta > 0,\gamma > 1}\sum _{i\in {{\mathcal{T}}}_{{\rm{bighit}}}}\log {p}_{\varDelta }({r}^{(i)}-{c}^{(i)})$$

This approach is conservative for our use of outlier detection as a null model. If the largest changes were not immediately followed by a decision, then delays would tend to be overestimated, causing the fitted outlier detection model to display longer-timescale dependencies that are typically associated with evidence accumulation strategies such as integration. Thus, the risk of falsely rejecting this null model would not increase.

For the largest change magnitude, miss trials predominantly reflected task disengagement rather than typical sensory/motor delays, and were therefore excluded when fitting the delay time distribution. According to a hidden Markov model, disengagement was the a posteriori most probable state for the majority of large-change miss trials (95.2% of large-change misses were during a disengaged state).

Fitting decision parameters. The decision parameters (sensory noise variance and decision threshold) were subsequently fit using the entire dataset, holding the delay time distribution fixed. Here in the general case, the decision and delay times cannot be directly observed, and were marginalized out as latent variables. The decision parameters were chosen to maximize the log marginal likelihood of the observed response data:

$$\mathop{\max }\limits_{{\sigma }^{2},b}\sum _{i\in {{\mathcal{T}}}_{{\rm{nonmiss}}}}\log {p}_{R}({r}^{(i)})+\sum _{i\in {{\mathcal{T}}}_{{\rm{miss}}}}\log (1-{F}_{R}({t}_{{n}^{(i)}}))$$

For hit and early lick trials (denoted ${{\mathcal{T}}}_{{\rm{nonmiss}}}$), the likelihood is given by the marginal probability density of a response at the observed movement onset time. For miss trials (denoted ${{\mathcal{T}}}_{{\rm{miss}}}$), the response time is treated as right-censored; its precise value is unknown, but is known to exceed the last time point in the trial. The likelihood for miss trials is thus given by the marginal probability mass lying beyond this point.

Sampling. To statistically compare mouse behaviour to the outlier detection null model, we sampled 10,000 synthetic datasets from the model fitted for each mouse. For every quantity of interest, the value computed from the real data was compared to values computed from each synthetic dataset, comprising 10,000 samples from the null distribution. Synthetic datasets were generated for each mouse as follows.

Each trial used the same change point and stimulus amplitudes presented in the real data. The real stimulus ended after the lick on trials where mice responded, leaving unknown future values that would have been presented had a lick not occurred. Such missing stimulus values were filled in by sampling from the same distribution used to produce the original stimuli (independently for each synthetic dataset).

Given the stimulus, a decision time and delay time were sampled from the distributions p_D and p_Δ described above. The sum of these quantities yielded a synthetic response time, representing movement onset.

To generate synthetic lick times, we assumed that the additional delay between movement onset and licking was i.i.d. across trials. We therefore sampled with replacement from the measured movement-to-lick delays in the real data. Synthetic lick times were obtained by adding sampled movement-to-lick delays to synthetic movement onset times.

Synthetic lick times were used to determine trial outcomes (hit, early lick, miss). Each trial was classified as a: hit if the lick occurred during the change period; early lick if the lick occurred before the change point; or miss if no lick occurred before the end of the change period.

Effect of magnitude and timing of TF pulses on probability of early licks

For analyses of the effect of TF pulses on probability of early licks we used the training data of the same 15 mice used for Neuropixels recordings. Here we only used sessions where mice reached robust proficiency of the task and were at the final training protocol (mean of 77.5 sessions per mouse). Note that here the time of lick onset was measured from the registration of lick by the spout as opposed to the videography analysis on Neuropixels recording sessions elsewhere in the manuscript. We used only trials where early licks happened at least 2 s after the baseline onset to decrease the influence of impulsive licks on results.

To empirically validate that mice use multiple pulses of sensory evidence to influence their decision to lick during the baseline period, we analysed how early lick probability is influenced by magnitudes and timing of preceding TF pulses. First, we tested whether the deviation of a single TF pulse relative to the mean baseline 1 Hz makes mice correspondingly more or less likely to make an early lick within the subsequent 0.2–1.0 s. For that we separated TF pulses by magnitude (in octaves) into 15 bins such that each bin contained approximately equal number of TF pulses. To calculate the conditional probability of early lick at a certain time after a TF pulse of a given magnitude, we found instances of such events (pulled across all sessions with robust performance for each mouse) and divided them by the total amount of early licks (Extended Data Fig. 6c). To calculate an overall influence of a TF pulse on early lick probability, we summed conditional probabilities within a [−1, −0.2] s window relative to early lick onset (Extended Data Fig. 6d):

$$P({\rm{L}}| {\rm{TF}})=\mathop{\sum }\limits_{t=-1}^{-0.2}P({\rm{L}}| {\rm{TF}}(t))$$

which can also be written as: P(L|TF) = P₀ + ΔP(L|TF), where

$${P}_{0}=\mathop{\sum }\limits_{t=-1}^{-0.2}P\left({\rm{L}}| {\rm{T}}{\rm{F}}\left(t\right)=1\,{\rm{Hz}}\right)$$

And can be thought as a chance level of making a lick without a deviation of stimulus TF from the mean baseline TF value.

The empirical effect of two TF pulses on lick probability was calculated from behavioural data in a similar way. To compare the measured effect of two TF pulses with their expected effect if they influenced the lick probability independently, we calculated their cumulative independent effect on early lick probability based on empirically measured effect of a single TF pulse on early lick probability. The independent effect of two TF pulses with a delay of Δt s between them can then be written as follows:

$${P}_{{\rm{Ind}}}({\rm{L}}| {{\rm{T}}{\rm{F}}}_{1},{{\rm{T}}{\rm{F}}}_{2})={P}_{0}+\Delta P({\rm{L}}| {\Delta {\rm{T}}{\rm{F}}}_{1})+\Delta P({\rm{L}}| {\Delta {\rm{T}}{\rm{F}}}_{2})-\Delta P({\rm{L}}| {\Delta {\rm{T}}{\rm{F}}}_{1})\times \Delta P({\rm{L}}| {\Delta {\rm{T}}{\rm{F}}}_{2})$$

where:

$$\Delta P\left({\rm{L}}| {\Delta {\rm{T}}{\rm{F}}}_{1}\right)=\mathop{\sum }\limits_{t=-1-\Delta t}^{-0.2-\Delta t}P\left({\rm{L}}| {\rm{T}}{\rm{F}}\left(t\right)\right)-{P}_{0}$$

$$\Delta P\left({\rm{L}}| {\Delta {\rm{T}}{\rm{F}}}_{2}\right)=\mathop{\sum }\limits_{t=-1}^{-0.2}P\left({\rm{L}}| {\rm{T}}{\rm{F}}\left(t\right)\right)-{P}_{0}$$

A deviation of lick probability after two TF pulses from the probability predicted by the independent effect of two TF pulses would indicate an interactive effect between pulses, which should be expected if mice utilize integration of sensory evidence. To measure the relative difference between the behavioural result and the expected independent effect of two fast TF pulses (Fig. 3d and Extended Data Fig. 6i), we calculated:

$$I=\frac{P({\rm{L}}| {{\rm{T}}{\rm{F}}}_{fast},{{\rm{T}}{\rm{F}}}_{fast})-{P}_{{\rm{Ind}}}({\rm{L}}| {{\rm{T}}{\rm{F}}}_{fast},{{\rm{T}}{\rm{F}}}_{fast})}{{P}_{{\rm{Ind}}}({\rm{L}}| {{\rm{T}}{\rm{F}}}_{fast},{{\rm{T}}{\rm{F}}}_{fast})}$$

When applying this analysis to the outlier detection agent data, we used data only from trials that resulted in early licks, meaning that the model made a decision to initiate a lick during the baseline period and before the TF change epoch. For outlier detection agent model that was fitted to a particular mouse data, we sampled the same number of early lick trials across 4,000 synthetic datasets (see section above) as there were present across all behavioural sessions of that mouse. The data was then pulled across all models corresponding to different mice and analysis steps were applied to the combined dataset as described above for the mice data. This procedure was repeated 4,000 times to estimate non-parametric 95% confidence intervals of results from the outlier detection agent.

Electrophysiological recordings

Prior to acute electrophysiological recordings, we habituated mice to the electrophysiological recording setup for 2–7 days (depending on the performance of the mouse in the electrophysiological recording setup), to allow mice to perform optimally during electrophysiological recording sessions.

Surgery

Once mice were habituated to the recording setup, we implanted a recording chamber with one or two 3 mm craniotomies inside, together with a stainless-steel grounding wire in the contralateral hemisphere, under 1.5% isoflurane together with administration of meloxicam (5 mg kg⁻¹) and dexamethasone (2–3 mg kg⁻¹). During surgery a kapton disk (Laser Micromachining Limited) was placed on top of the dura inside each craniotomy. The disk had 19 holes with 0.5 mm diameter, arranged in a honeycomb shape, for keeping track of probe insertions. The craniotomy and disk were covered with DuraGel (Cambridge NeuroTech) to protect the brain. A 1–2 mm tall plastic enclosure was then positioned around craniotomies and sealed around the edges with bone cement. Finally, we covered the plastic enclosure with a removable plastic cover, to create a rigid physical barrier over the DuraGel sealed craniotomy, to provide robust protection of the recording preparation between recording sessions. The mice were allowed to recover for 24 h before the first recording session took place.

Recordings

Electrophysiological data collection was done using Neuropixels 1.0 probes (IMEC, Belgium) and collected with a PXI based system (National Instruments), and saved using SpikeGLX (https://github.com/billkarsh/SpikeGLX). For trained mice, we recorded up to 13 sessions per mouse (167 probe insertion from 114 sessions total (15 mice)). For untrained mice, we recorded up to 9 sessions per mouse (89 probe insertions from 45 sessions total (6 mice)). Probes were dipped in CM-DiI (Sigma-Aldrich) prior to insertion. In each session, we inserted up to 2 probes at a time. The probes were always inserted at the same angle within the coronal plane (10° and −15° relative to the vertical axis) to aid subsequent histological probe tract tracing.

At the beginning of each session, we removed the plastic lid above the recording chamber exposing the DuraGel covered craniotomy, and inserted the probe(s) through the DuraGel using microcontrollers (Sensapex) at 5–10 μm s⁻¹. The probe(s) was allowed to settle for 20 min, to increase stability throughout the recording session. At the end of the session probes were removed (at 15 μm s⁻¹) and the plastic cover over the recording chamber was reattached for protection of recording preparation.

The setup for presenting stimuli and monitoring behaviour were identical to the setups in which mice had been trained (see ‘Behavioural task’).

Pre-processing and spike sorting of electrophysiological data

Electrophysiological data was first filtered using CatGT (https://billkarsh.github.io/SpikeGLX/#catgt) with modified form of common average referencing (-dlbdmx flag).

Spike sorting. We spike-sorted electrophysiological data from each probe in each session using KiloSort2.0⁶⁵ (https://github.com/MouseLand/Kilosort). For initial selection of units undergoing further curation, we only selected units designated as ‘good’ (based on cross-correlogram contamination) by KiloSort2.0.

Quality checks. For our electrophysiological recordings of trained mice, we manually inspected and curated, in Phy2.0 (https://github.com/kwikteam/phy), every unit which KiloSort2.0 had designated as ‘good’. For our recordings in trained mice this left 44,288 units to be manually inspected and curated, and 15,406 units were kept for analysis after manual curation. Based on the manual curation data from trained mouse recordings (see ‘Manual curation of spike-sorted units from trained mice’), we established a series of heuristics for creating automatic curation of units (see ‘Automatic curation of spike-sorted units from untrained mice’) and used these for recordings from untrained mice.

Manual curation of spike-sorted units from trained mice. We manually inspected and curated all units which KiloSort2.0 had designated as good, based on cross-correlogram contamination. In Phy2.0, we first inspected and merged units that clearly belonged to the same cluster, but had been split by KiloSort2.0, or split the noise from signal in units with clearly separatable noise contamination. We then designated each unit into one of five categories:
1. (1)
  
  Perfect, or almost perfect, with no/very minimal noise, drifting, cutting in/out for the full duration of recording.
2. (2)
  
  Usable and good signal with some noise that cannot be extracted that lasts for the full duration of the recording.
3. (3)
  
  Some drift, but possible physiological change in signal. Clear signal for most of duration of the recording.
4. (4)
  
  Drifting/sudden loss, but otherwise usable/close to perfect. Clear signal for over 50% of the duration of the recording but requires only using a subset of the session.
5. (5)
  
  Noise/useless. Spike shape is not physiological.
Our goal was to remove from analyses units that had large contamination with multi-unit activity, were not recorded throughout the full duration of a session, or were a result of artifacts in recorded signals. We therefore used units designated as category 1–3 above for all further analysis from trained mice.

Automatic curation of spike-sorted units from untrained mice. We next used the manual designations of units to establish a set of criteria for automatic detection of units we would include with manual curation. Based on the manual curation data above we established the following 7 criteria for considering a unit good for analysis:

Firing rate criteria:
1. (1)
  
  Mean firing rate must be above 0.5 Hz.
2. (2)
  
  Rolling 20-min average firing rate cannot drop below 30% (that is, 70% drop from mean) of its mean firing rate.
3. (3)
  
  Rolling 10-min average firing rate cannot drop below 20% (that is, 80% drop from mean) of its mean firing rate.
4. (4)
  
  Rolling 5-min average firing rate cannot drop below 10% (that is, 90% drop from mean) of its mean firing rate.
5. (5)
  
  Inter-spike interval (ISI) violations. Absolute refractory period needs to have <20% estimated contamination rate from other neurons (this is what Kilosort2.0 calls ‘good’).
6. (6)
  
  If there are some spikes in the refractory period, the ISI peak in the first 5 ms cannot be within the first 2 ms.
7. (7)
  
  ISI histograms cannot have sudden large spikes in their shape (that is, peak of ISI cannot be 4 times larger than the second highest peak—that is usually its immediate neighbour).
These criteria selected approximately 90% of units we would have designated with categories 1 (perfect, or almost perfect) or 2 (usable and good signal with some noise) with manual curation, and excluded approximately 85% we would have designated as 4 (drifting/sudden loss) or 5 (noise) with manual curation.

This automatic selection of units was used to select units for analysis from untrained mice recordings and yielded 6,215 units out of 20,292 ‘good’ KiloSort2.0 units.

Clock-drift correction. A shared 1 Hz square wave signal was recorded on the clock of each headstage and National Instruments (NI) acquisition card using a SYNC option in SpikeGLX. Clock drift between spike times from different probes and behavioural events extracted from NI acquisition card recording was corrected post-hoc via TPrime (https://billkarsh.github.io/SpikeGLX/#tprime) using the shared square wave signal.

Videography

Acquisition

High-speed videography of front (100 frames s⁻¹, 640 × 512 pixels) and side view (50 frames s⁻¹, 976 × 1,024 pixels) of the mouse face was acquired using two Chameleon3 cameras (CM3-U3-13Y3M-CS, FLIR) with infrared illumination. The videos were acquired in an 8 bit greyscale format. Cameras were configured to send a TTL signal to the National Instruments PXIe board at the start of exposure of every acquired frame. These TTL signals were used to align frame times to the time of behavioural events and spike times.

Pupil size

In order to estimate the pupil size, we trained DeepLabCut⁶⁶ to track the pupil size and position using videos acquired with the side camera. The model was trained to track 12 points surrounding the mouse pupil. In order to assess the model performance, after the training the model was tested on videos from sessions not used for training. Pupil size was estimated as an area of an ellipsoidal best fit to the tracked 12 points surrounding the pupil.

Motion energy

For calculation of motion energy, we primarily used videos acquired with the front camera to access a finer temporal resolution (with the exception of 2 sessions where for technical reasons we used a lower jaw ROI from side camera video). To estimate motion onset times, we used ROI centred around the mouse’s face, though nearly identical results were obtained with lower jaw or whisker pad ROIs from the side camera (data not shown). Motion energy was defined as a square root of the sum of squared frame-to-frame pixel value differences, divided by the number of pixels within the ROI.

Movement onset time estimation

In order to find the onset times of orofacial movements, we wanted to estimate the typical noise level of the motion energy signal and find the time points where the signal significantly deviated from the noise-band level. As a first step, we calculated the distribution of motion energy values in a 2-s window centred around the lick registration times. We next fitted a mixture of Gaussian distributions with the goal to capture both contribution of the variance of motion energy values during the lick as well as due to noise. The mixture of three Gaussian distributions worked well to fit the data across all sessions and mice. The threshold for the presence of movement was defined as the mean plus two standard deviations of the Gaussian with the lowest value of the mean from the Gaussian mixture.

Finally, to find the time of motion onset time, we looked backwards in time from the time of lick registration by the piezo signal. The time point preceding the first instance of motion energy going below the threshold value defined above was considered the onset time of the orofacial movement.

Histology

For histological identification of the location of the recording probes and allocation of unit location in the mouse brain, we followed a protocol similar to ref. ⁶⁷.

Serial 2-photon tomography for Neuropixels probe tract tracing

Following a terminal administration of pentobarbital, mice were perfused with a phosphate buffer solution (PBS) followed by 4% paraformaldehyde (PFA) solution. We post-fixed the brain in the 4% paraformaldehyde for a minimum of 24 h at approximately 5 C. Following fixation, brains were moved to PBS for a minimum of 12 h prior to imaging. For imaging, brains were embedded in 5% agarose gel and mounted onto a vibratome cutting stage under the microscope objective. The brains were imaged using serial section two-photon microscopy⁶⁸. The microscope was controlled with ScanImage Basic (Vidrio Technologies), and custom software (BakingTray (https://github.com/SainsburyWellcomeCentre/BakingTray)). Images were stitched into a full 3D rendering of the brain using custom software (StitchIt (https://github.com/SainsburyWellcomeCentre/StitchIt)). We imaged the entire brain (from the olfactory bulb to the beginning of the spinal cord) with a resolution of x: approximately 2 μm, y: approximately 2 μm, z: 20 μm, with a 920 nm two-photon laser (100–150 mW power at sample). We sliced the brain in 40-μm sections, and imaged 2 z-planes (around 25 μm and around 45 μm from the tissue surface) into the remaining tissue following each 40-μm section. Two PMTs, one for capturing green (bandpass filter ET525/50 m) and red (bandpass filter ET570lp) fluorescence acquired the 2 channels of data subsequently used for analysis.

Neuropixels probe tract alignment to the Allen Common Coordinate Framework atlas and estimation of unit location

Prior to image processing, we downsampled microscopy images to 10-μm voxels and registered the brain to the standardized Allen Common Coordinate Framework (Allen CCF⁶⁹) using custom software (BrainRegister (https://github.com/stevenjwest/brainregister)). We then manually traced each neuropixels probe tract through the brain in 3D using custom software (Lasagna (https://github.com/SainsburyWellcomeCentre/lasagna)). Finally, we assessed the overall firing rates and LFP spectra of individual Neuropixels channels and compared it to atlas positions. Where needed, we manually adjusted the scaling of brain regions along the probe track to align responses on channels with features associated with anatomical locations using custom software (Ephys alignment tool (https://github.com/int-brain-lab/iblapps/tree/master/atlaselectrophysiology)⁷⁰). Unit location was estimated from the location of the channel that had the largest absolute peak value of the mean waveform. For all analyses, we combined units across all subdivisions of a brain region (layers of cerebral cortex, dorsal and ventral divisions as ACAd and ACAv and in some cases functionally similar brain regions—see Supplementary Tables 1 and 2).

Neural data analysis

Only brain regions with at least 40 units were analysed. Analyses specific to TF-responsive units were done only for brain regions with ≥10 such units. No further sample size calculations were performed. Manual curation of units’ quality and stability was done without the knowledge of brain regions from which recordings were made. The subsequent analyses pipeline was applied in the same manner to data from all applicable brain regions, but the custom nature of analyses prevented investigators to remain blind to the identity of brain regions or dataset type (trained versus naive mice).

GLM of neural activity

Model. We binned neural activity in 50-ms bins (matching the duration of each TF pulse) aligned to trial start. We then fitted a Poisson generalized linear model to predict trial-to-trial neural activity as a function of a set of temporally unfolded task-related predictors that were present during a trial. Each predictor was extended temporally prior and/or post the timing of the predictor in 50-ms discretized steps (matching neural activity binning), with an independent weight estimated for each time step around the predictor. We predicted neural activity using 19 task-related predictors:

(1) TF fluctuations during baseline period (kernel length: 0–1.5 s); (2) Trial start (0–1 s); (3) Time since baseline start (from 1 s from trial start to change onset); (4–9) Six change onsets (a separate predictor for each change size (0–2 s)); (10) Lick preparation (−1.25–0 s prior to lick); (11) Lick execution (0–0.5 s post lick); (12) Air-puff (0–0.25 s); (13) Reward (0–0.4); (14) Abort (−1.25–0.25); (15) Phase of grating for upwards drift (12 phase bins from 0–360°); (16) Phase of grating for downwards drift (12 phase bins from 0–360°); (17) Video motion energy (−0.05–0.8 s); (18) Running wheel movement (−0.05–0.8 s); (19) Pupil diameter (−0.75–0.75 s).

We fit the model with L2 (ridge) regularization, optimized with cyclical coordinate descent as implemented in GLMnet⁷¹ (α = 0). We trained a model for each neuron on 90% of the data, and cross-validated on 10% of the data, and iterated the predictions over a tenfold cross-validation. Within the training dataset we tuned the L2 regularization term using tenfold cross-validation.

Identification of units encoding TF, lick preparatory activity and/or lick execution activity. To identify which cells significantly responded to a predictor of interest (that is TF fluctuations during baseline, lick preparation epoch, or lick execution epoch), we first re-fitted reduced models similar to the full model on 90% of the data, with 10-fold cross-validation, except we removed a predictor(s) of interest: (1) For identification of TF-responsive units, we estimated a model where we removed the predictor estimating the responses to TF fluctuations during baseline. (2) For identification of units with lick preparation activity, we estimated a model where we removed the predictor estimating the activity leading up to a lick. (3) For identification of units responding to lick execution, we estimated a model where we removed the predictor estimating activity during lick execution, the predictor estimating activity modulation by motion energy captured by videography, and the predictor estimating activity modulation by running wheel movement.

For each 10% test set, for each neuron we then calculated the mean actual peri-event time histogram (PETH) as well as the mean predicted PETH of both the full model and the reduced model for the following types of events: (1) −0.15 to 0.75 s around fast and slow TF pulses (that is, TF values 0.5 s.d. from the mean TF during baseline); (2) −1.5 to 0 s prior to early lick onsets; and (3) 0 to 0.4 s post lick onset.

A unit was considered significantly encoding TF pulses during the baseline period if two criteria were satisfied: (1) The mean Pearson’s correlation prediction of the full model (across k-folds) from the combined mean fast and slow TF pulse response (that is, mean fast TF pulse and mean slow TF pulse responses subtracted from each other) was >0.2; and (2) if the cross-validated prediction of the TF response after subtracting the predicted TF response of the reduced model with no TF fluctuation predictor—that is, residual prediction—was significant (P < 0.01 (t-test), n = 10 independent cross-validations). A unit was considered significantly encoding lick preparation if (1) the mean Pearson’s correlation prediction of the full model (across k-folds) of the mean activity leading up to a lick (−1.25 to 0 s) was >0.2; and (2) if the cross-validated prediction of the mean activity after subtracting the predicted mean activity of the reduced model with no lick preparation kernel—that is, residual prediction—was significant (P < 0.01 (t-test), n = 10 independent cross-validations). Finally, a unit was considered significantly encoding lick execution if (1) the mean Pearson’s correlation prediction of the full model (across k-folds) of the mean activity following a lick (0 to 0.25 s) was >0.2; and (2) if the cross-validated prediction of the mean activity after subtracting the predicted mean activity of the reduced model with no lick preparation kernel—that is, residual prediction—was significant (P < 0.01 (t-test), n = 10 independent cross-validations).

Focality index. To assess how distributed TF encoding was across brain areas, before and after learning, we computed a focality index (F) (similar to Steinmetz et al.³) of the TF encoding:

$$F=\frac{\sum ({p}_{a}^{2})}{{(\sum {p}_{a})}^{2}}$$

where p_a is the proportion of neurons in an area that is encoding stimulus TF during the baseline period. If all TF encoding neurons were confined to a single area, this measure would take on the value of 1. If encoding was perfectly distributed across all areas recorded this measure would take on the value 1/N_areas. In order to compare between untrained and trained mice, we identified the common areas which had more than 40 units recorded in both trained and untrained mice. This left N = 24 areas from which to estimate the focality index. We estimated 95% confidence intervals and P values by bootstrapping the neurons included in the estimation 10,000 times with replacement.

Peak time and width of GLM estimated TF kernels for TF-responsive neurons. To investigate the peak time and width of the GLM estimated TF kernel for assessing how sustained responses to TF fluctuations were based on GLM weights, we first identified the absolute peak value of the TF kernel; because the GLM was based on 50 ms binning of spike counts, peak times for the GLM TF kernel was in 50 ms resolution. In cases where the absolute peak position within 1 s was a negative weight, we flipped the kernel in order to calculate the width. We then estimated the full width at half maximum (FWHM) of each TF kernel around its peak using findpeaks in MATLAB. For each area, we calculated the median peak time and median FWHM across all TF-responsive units.

Ramping differences in GLM change kernels. To test how neurons accumulated evidence when they were presented with a rewarding sustained change in stimulus speed, we tested how the slope of the visual evoked ramping activity following a change onset was dependent on the amount of evidence (change size) being presented. To isolate the visual component of the activity following change onset, we used the GLM kernel which fits the activity following change onset until change offset, while linearly taking into account other variables which may contribute to activity such as pupil size, preparatory activity and movement-related activity (see Model).

We estimated the mean change kernel for each change size for TF-responsive and non-TF-responsive units separately for each area. In cases where responses to fast TF pulses were negative, we flipped the change kernel so every unit had responses aligned to positive fast TF pulses—this allowed the mean to capture the visual evidence activity ramp irrespective of sign. We then identified the time point for each change size where the change kernel reached 50% of its maximum weight (To control for noise fluctuations in kernel weights, we approximated the 50th percentile crossing by taking the mean time point of the 33.33rd percentile, 50th percentile and 66.66th percentile crossing). We then calculated the degree to which activity ramping time scaled with change size, by regressing the 50th percentile crossing against change size. We estimated the non-parametric 95% confidence intervals and P values of the relationship between change size and 50th percentile crossing (that is, ramping time/change size) by bootstrapping with replacement (10,000 times) the neurons went into the mean change kernels, and then estimating the slope of the regression for each bootstrapped mean change kernels.

Propagation and widening of TF pulse evoked activity

Identification of TF pulse outlier events. Fast TF pulse was defined as TF fluctuations larger than 1 s.d. of baseline TF fluctuations (in log₂ scale) above the mean TF value (TF > 1.19 Hz). Similarly, slow TF pulse was defined as TF fluctuations below 0.84 Hz.

For calculation of average response to TF outlier events, we considered only TF outlier events satisfying the following criteria:
1. (1)
  
  Later than 1 s from the baseline onset.
2. (2)
  
  Earlier than 2 s + post pulse analysis window from the motion onset time on early lick or abort trials.
3. (3)
  
  Excluding the change period plus a post pulse analysis window.
The aim of these criteria was to exclude the influence of baseline onset, movement, or preparatory activity on the response to TF pulses.

Estimation of peak time and width of TF pulse evoked activity. For each unit defined as TF-responsive by the GLM analysis described above, we calculated a mean response to a fast pulse using outlier events that occurred during the baseline period and satisfied the criteria outlined above. Additionally, we calculated a mean response to TF pulses within [−0.5, 0.5] s.d. of the baseline TF fluctuations. The goal of this procedure was to capture continuous ramps of activity that some units exhibited and exclude their influence on the shape of response to a TF pulse. We applied the subtraction of this baseline response for all TF pulse response analysis unless explicitly stated.

Next, for the baseline subtracted mean response to a fast TF pulse, we calculated its peak time, as the time of the largest absolute change in firing rate within 1 s from the pulse onset, and a corresponding half-peak width.

Integration of multiple TF pulses. Because the noise in TF fluctuations is random, by chance there are occurrences of two fast pulses separated by a certain delay. To study the integration of TF pulses, we found such instances of events where two fast pulses occurred at a given delay between the offset of the first and the onset of the second, additionally also satisfying the exclusion criteria outlined above. The mean response aligned to such events was considered a response to a sequence of two fast pulses.

For computing the mean response across all TF-responsive units within a brain region, in order to avoid averaging across responses with different signs, we flipped the sign of response for units that showed decreases in activity after a single fast pulse. For computing a z-score of response, the mean and s.d. were estimated from 0.5 s preceding the first pulse onset.

Facilitation by the second fast pulse. First, we measured an average of z-scored responses across the population of TF-responsive units within a brain region to a single fast TF pulse. We then computed the peak value of that response (r_1fast), and a corresponding peak time. To find the size of response to a sequence of two fast pulses (r_2fast), we found a time point at the same delay from the onset of the second fast pulse as the peak time of response to a single fast pulse and found a peak value of response within 100 ms centred around that time point. The relative facilitation to a sequence of fast pulses was defined as $\varDelta \,=\,\frac{{r}_{2{\rm{fast}}}\,-\,{r}_{1{\rm{fast}}}}{{r}_{1{\rm{fast}}}}$.

To determine the confidence intervals for the results of this analysis, we bootstrapped with replacement (2,000 times) across TF-responsive neurons and repeated the analysis described above for each sample of neurons. Shaded regions indicate 2.5 and 97.5 percentiles of the resulting distribution.

Preparatory activity before the lick onset

To study change-aligned (Fig. 3) or hit lick-aligned (Fig. 5) activity, we computed z-score of mean PETH for each unit. z-Scoring was done using the mean and s.d. estimated from activity during 2 s before the change onset.

For analysis shown on Fig. 5, for each brain region the fraction of significantly active units within a group (that is, TF-responsive) was measured by calculating at every time point a fraction of units with the absolute value of z-score larger than the significance threshold of 2.576 (corresponding to P < 0.01). Additionally, we subtracted the ‘baseline’ level of activity calculated within [−2, −1.8] s before hit-lick onset, which for a few brain regions was larger than chance level likely due to non-normal distribution of firing rates or a small number of events used for estimation of the mean and s.d. The confidence intervals were estimated by bootstrapping with replacement (5,000 times) across TF-responsive (or TF non-responsive) neurons and repeating the estimation of fraction of significantly active neurons for each sample of neurons.

The latency of activation of TF-responsive or TF non-responsive populations was defined as the earliest time point following which within a 100-ms window for at least 80 ms: (1) the lower 95% confidence interval of fraction of active units was above zero; and (2) the mean fraction of active units was above 0.1.

The latency of significant difference in activation between TF-responsive and TF non-responsive populations was estimated as the first time point where within a 100-ms window for at least 80 ms the confidence intervals of the difference in activation were above zero.

The latency of significant difference in activation across all units in each brain region (Extended Data Fig. 9a) was estimated as the first time point where within a 100-ms window: (1) the lower 95% confidence interval of fraction of active units was above zero; and (2) the mean fraction of active units was above 0.05.

Intrinsic timescales

We binned the neural activity into 50-ms bins (same binning was used in ref. ²³). We then calculated the temporal autocorrelation (20 lags = 1 s) of spike counts using Pearson’s correlation in the inter-trial intervals between −2.5 s to −0.5 s prior to trial onset for each neuron (in this period mice were seeing a grey screen, and trained mice had to remain stationary for at least 3 s for the trial to begin).

To determine the intrinsic timescale for each area, we fit an exponential decay function to the mean autocorrelation function of all the units recorded in the area. For single-neuron analysis of relationship between intrinsic timescales and TF width, we estimated the autocorrelation for each TF-responsive neuron separately. For areas or neurons with autocorrelation functions with non-monotonic decay, we fit the exponential decay from the part of the autocorrelation where monotonic decay was happening (in a subset of areas this would mean offsetting the fit 1–3 time bins). Finally, we calculated the τ (that is, the intrinsic timescale value) of the exponential decay (accounting for offset where necessary).

Population analysis

Similarity of TF pulse responses and lick preparation activity in TF-responsive populations

We assessed the similarity of TF responses and lick preparation activity across TF-responsive populations in each area by estimating the Pearson’s correlation of mean firing rates (within a 50-ms window around the mean activity peak time across neurons within the each area) following fast TF pulses (that is, >1 s.d. TF value) and their mean activity prior to early lick [−0.3 to 0 s] (after normalizing firing rates by subtracting baseline firing rates from both TF responses and lick preparation activity). We estimated the non-parametric 95% confidence intervals and P values by bootstrapping with replacement (10,000 times) the neurons going into the correlation.

Pre-processing steps

For all units located within a given brain region, but not necessary simultaneously recorded, we first computed the mean neural responses across a given trial type (for results shown on Fig. 6b–h: hit trials during weak TF changes (1.25 and 1.35 Hz) aligned to the lick onset times, [−2, 1.5] s time window). Only trials with hit-lick onset times larger or equal 0.4 s from change onset were used. Neurons from sessions with less than 10 trials of a given type were excluded from this analysis. Firing rates were calculated as spike counts averaged in 10 ms bins and smoothened by convolution with two-sided Gaussian with 30 ms s.d. The mean neural responses were combined into a firing rate matrix (but also see cross-validation section) with dimensions of Neurons × Time.

Neural data was pre-processed in the following way: first, to limit the dominant influence of high-firing units, we applied soft-normalization to each neuron’s firing rate, such that the neurons with strong responses had close to unity range of responses ${r}^{/}\,=\,\frac{r}{7+(\max (r)-\min (r))}$. The constant 7 was chosen as the roughly 20th percentile value of the firing rate range across all units. Second, the neural responses were mean-centred by subtracting the mean of each neuron’s activity across time and the mean activity across all neurons at every time point.

Definitions of movement and movement-null subspaces

We used the approach first utilized in ref. ³³. There, the authors formalized a method to find a linear mapping between low-dimensional representation of activity in PMd/M1 and the muscles EMG data, which defines a movement subspace. A null-space relative to that subspace forms an orthogonal set of dimensions which activity can occupy without directly affecting the movement execution. To extend this analysis on our data, we used combined recordings of orofacial motor and premotor nuclei (V, IRN, SPVI and SPVO) as a proxy for activity of orofacial muscles involved in execution of a lick. While recordings from GRN could have also been included into this group, we kept it separate to allow the population analysis to be applied to that region because (1) we had a large number of units recorded from that region alone; and (2) it was the only nucleus in medulla with above-chance number of TF-responsive units, warranting a separate analysis.

We considered a possible mapping onto the movement subspace for each brain region. Our rationale was the following: there exist several parallel neural pathways that can drive the activity of orofacial nuclei neurons–from primary motor cortex, basal ganglia, cerebellar or midbrain output regions^56,57,72. Thus, the modes of activity within these regions that map onto the movement subspace may have a causal role for the execution of licks. In general, however, these signals can also be caused by movement afference that is broadcasted globally^3,5,9 (Fig. 1k,l). It is impossible to differentiate between these two possibilities from our data alone and thus the existence of mapping of activity onto the movement subspace does not necessarily imply that the brain region is causally involved in execution of the lick. With that said, we did not find a good mapping onto a movement subspace for most of the early visual areas, olfactory regions and hippocampal input regions (Extended Data Fig. 10a), suggesting that existence of mapping onto the movement subspace is not possible across all brain regions.

The mapping onto movement subspace was defined as:

$$\widetilde{M}=W\widetilde{N}$$

(1)

where $\widetilde{M}$ and $\widetilde{N}$ are low-dimensional representations of activity (projections onto main the principal components, the latter found via svd Matlab function) of neurons within the orofacial nuclei group and the target brain region, respectively, and W is a linear mapping operator onto the movement subspace.

Before finding a linear mapping, we also zeroed the initial state across projections on principal components by subtracting from each projection the mean value within [−2, −1.5] s from lick onset. This step avoided the need for using intercept in the linear fit and simplified the visualization of projections on principal components and movement/movement-null dimensions. Linear mapping was found using only the time-period containing movement-related activity of orofacial nuclei [−0.1, 1.5] s around lick onset. This way we did not preclude the presence of preparatory activity on movement dimensions from the definition of the linear fit itself. A linear mapping to movement dimensions was found using linear regression with the Matlab function lsqnonlin.

Correspondingly, W_null was a null-space of W and was found using the Matlab function null. We used two top principal components of orofacial nuclei activity (which captured 61% of the total cross-validated variance; Extended Data Fig. 10a,b) and 4 top principal components of activity in a target brain region to find W and W_null operators (see Extended Data Fig. 10b–d). This choice resulted in both movement and movement-null subspaces being two-dimensional. We additionally ensured that norms of these operator are equal $\parallel {W}_{{\rm{null}}}\parallel =\parallel W\parallel $ in order to make the comparison between the movement and movement-null subspaces fair.

Since the definition of specific dimensions in movement-null subspace is to a degree arbitrary, we defined the first movement-null dimension by finding a rotation within the movement-null subspace that maximized the amount of variance captured by that dimension prior to lick onset. The second movement-null dimension was then simply orthogonal to the first dimension in movement-null subspace. This was used mainly to simplify visualization, with all subspace-related analyses done using both dimensions in each subspace.

The positive direction of movement dimensions was chosen such that the mean value of projection of orofacial nuclei activity within [−2, 0.5] s around lick onset was positive. The positive direction for movement-null dimensions was chosen such that the mean value of projection of activity within [−2, 0] s around lick onset was positive.

Subspace occupancy

Relative subspace occupancy at a moment of time t was defined as

$${O}_{{\rm{R}}}(t)=\frac{{E}_{{\rm{null}}}(t)-{E}_{{\rm{m}}}(t)}{{E}_{{\rm{null}}}(t)+{E}_{{\rm{m}}}(t)}$$

where E_null(t) and E_m(t) are Euclidean distances within movement-null and movement subspaces, measured between the neural state at the current moment of time t and the initial time point (the mean across 2 and 1.5 s before the lick onset). Values close to zero signify equal occupancy between subspaces and positive values indicate a preferential occupancy of the movement-null subspace. The peak-normalized occupancy (Extended Data Fig. 11a,b) was defined as $O(t)=\frac{E(t)}{\max (E)}$

Decomposition of projections onto contributions from TF-responsive and TF non-responsive units

We decomposed the projections on main principal components into a sum of contributions from TF-responsive and the TF non-responsive units. For that, we used the knowledge of identity of each unit as TF-responsive or TF non-responsive and wrote down the principal components U (from the singular value decomposition (SVD) of the firing rate matrix N = USV^T) as a sum of two parts as:

$$U={U}_{{\rm{TF}}}\left(\begin{array}{c}{U}_{i},i\in {\rm{TF}}\,{\rm{unit}}\\ 0,i\in {\rm{non}}\,{\rm{TF}}\,{\rm{units}}\end{array}\right)+{U}_{{\rm{non}}}\,{\rm{TF}}\left(\begin{array}{c}0,i\in {\rm{TF}}\,{\rm{unit}}\\ {U}_{i},i\in {\rm{non}}\,{\rm{TF}}\,{\rm{units}}\end{array}\right)$$

(2)

where U_i is a loading of the ith unit.

With that, projections on principal components can be written as:

$$\widetilde{N}={U}^{T}N={U}_{{\rm{TF}}}^{T}N+{U}_{{\rm{non}}}\,{{\rm{TF}}}^{T}N={\widetilde{N}}_{{\rm{TF}}}+{\widetilde{N}}_{{\rm{non}}}\,{\rm{TF}}$$

(3)

Substituting equation (3) into equation (1) gives projections onto movement dimensions as:

$${\widetilde{N}}^{{\rm{mov}}}=W\widetilde{N}={\widetilde{N}}_{{\rm{TF}}}^{{\rm{mov}}}+{\widetilde{N}}_{{\rm{non}}}\,{{\rm{TF}}}^{{\rm{mov}}}$$

and, correspondingly, projections on movement-null dimensions are written as:

$${\widetilde{N}}^{{\rm{null}}}={W}_{{\rm{null}}}\widetilde{N}={\widetilde{N}}_{{\rm{TF}}}^{{\rm{null}}}+{\widetilde{N}}_{{\rm{non}}}\,{{\rm{TF}}}^{{\rm{null}}}$$

The relative contribution of TF-responsive units within movement and movement-null subspaces at the moment of time t was then defined as following:

$$\begin{array}{l}{w}_{{\rm{mov}}}(t)=\frac{{\widetilde{N}}_{{\rm{TF}}}^{{\rm{mov}}}(t)\times {\widetilde{N}}^{{\rm{mov}}}(t)}{{\parallel {\widetilde{N}}^{{\rm{mov}}}(t)\parallel }^{2}}\times \frac{{\widetilde{N}}^{{\rm{mov}}}(t)}{| {\widetilde{N}}^{{\rm{mov}}}(t)| }\\ {w}_{{\rm{null}}}(t)=\frac{{\widetilde{N}}_{{\rm{TF}}}^{{\rm{null}}}(t)\times {\widetilde{N}}^{{\rm{null}}}(t)}{{\parallel {\widetilde{N}}^{{\rm{null}}}(t)\parallel }^{2}}\times \frac{{\widetilde{N}}^{{\rm{null}}}(t)}{| {\widetilde{N}}^{{\rm{null}}}(t)| }\end{array}$$

where the second multiplicative term ensures that the sign of contribution is relative to the defined positive direction (see above) of dimensions within each subspace.

In order to test whether the contribution of TF-responsive units is larger than what is expected from a uniform contribution of the full population, we repeatedly randomly selected (2,000 times) the same number of units as there were TF-responsive ones from the whole population and computed their contribution to projections on movement and movement-null dimensions as described above.

In addition to the analysis described above, we have also checked whether the above-chance contribution of TF-responsive units is a consequence of their level of activity, despite the normalization method that we used, or does it reflect a better correspondence of their activity to the population modes of activity within the movement-null subspace. For that we looked at the distribution of loadings along the first movement-null dimension–that captured the majority of preparatory activity there. We found that the majority of brain regions where TF-responsive units had above-chance contribution to the preparatory activity also had larger absolute values of loadings along that dimension than the rest of the population (Extended Data Fig. 11d,e).

Cross-validation

Since our analyses were focused on characterizing the mean neural responses, the cross-validation procedure that we used was designed to test the stability of the mean neural responses and their corresponding low-dimensional representations across trials. For that, we split trials into two randomly assigned and equally sized groups (fit and test trials) and calculated the mean neural response per unit across each group of trials. We next combined firing rates of neurons from the same brain region(s) (but not necessarily simultaneously recorded) into a joint matrix. After applying the pre-processing steps outlined above, we had two firing rate matrices from fit and test trials.

For cross-validated PCA (Extended Data Fig. 10a), we applied SVD on the first (fit) matrix and measured how well the remaining (test) matrix is predicted by the reconstruction from SVD components found from the first matrix. Similarly, the projections of activity on main principal components (Extended Data Fig. 13) were done using the test data, projected onto principal components found from the fit data.

For further analyses utilizing movement and movement-null subspaces, we applied SVD separately on each matrix and found their projections on first four main principal components. We then used low-dimensional representation of fit trials data to find linear mapping W and W_null onto the movement and moment-null subspaces. Finally, we applied W and W_null found from the fit data to the low-dimensional representation of the test data. This procedure was repeated 2,000 times, the 95% confidence intervals shown in Fig. 6 illustrate the 2.5 and 97.5 percentiles across projections of the test data. Because the sign of projection is arbitrary defined, we additionally applied a potential flipping of the sign of eigenvectors from each draw based on which direction had better alignment with the eigenvectors computed from the full firing rate matrix without the split into fit and test trials.

Responses to TF pulses

For each brain region, we constructed a firing rate matrix of all units responses to a fast TF pulse (or concatenating in time responses of each unit to different types of TF pulses for analysis shown in Fig. 6i,k–m), and used the same pre-processing steps as described above. The projections onto the movement and movement-null dimensions were done using loadings found from the analysis of hit licks activity described above (using the full firing rate matrix of hit-lick responses without the split into fit and test trials). Cross-validation of consistency of projections was done by randomly selecting half of TF outlier events, computing the mean firing rate across those events for each unit, applying the steps above to find the projections, and repeating this procedure 2,000 times. For analyses where different brain regions were combined into a common group, all units from those brain regions were combined into a joined firing rate matrix and the steps described above were applied.

Alignment of fast TF pulse response with a given dimension in movement or movement-null subspace was calculated as a cosine of an angle between the projection onto a target dimension and a 4-dimensional vector of TF pulse response (2 movement and 2 movement-null dimensions) at a time of the maximum Euclidean distance from the initial state across 4 dimensions within a 0.75-s window from the pulse onset. Similarly, for calculating the scaling of responses to different TF pulses along the first movement-null dimension, we found the sizes of projections at times of maximal Euclidean distance from the initial state within a 0.75-s window from the first TF pulse onset.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
[ad_2]

Source link
September 11, 2024