Dynamical structure of cortical taste responses revealed by precisely-timed optogenetic perturbation

The purpose of perception is to drive action - during tasting, for instance, every stimulus necessarily drives either swallowing or ejecting (“gapes”). Taste responses in the rodent primary gustatory cortex (GC) span this sensorimotor divide, progressing through a series of firing rate epochs that culminate in the emergence of palatability-related firing. Population analyses reveal palatability-related activity to appear through sudden and coherent ensemble transitions that, despite occuring between 0.5s and 1.5s on individual trials, reliably precede gaping onset by 200-300ms. Here, we tested whether these transitions actually drive gaping, delivering 0.5s perturbations of GC at various points in tasting trials. Perturbations significantly delayed gaping, but only when they arrived before palatability coding. Thus, perturbation had no impact on trials in which the transition had already occurred, but the identical perturbation delayed gaping on trials in which it hadn’t. Our results suggest a distributed attractor network model of taste processing, and a dynamical role for cortex in driving motor behavior.


Introduction
One of the primary roles of sensory processing is to drive action, in order that the source of sensory information can be either acquired or avoided (Prinz (1997), Wolpert and Kawato (1998), Wolpert and Ghahramani (2000)). Implicit in this statement is a suggestion that sensory and motor processing are tightly coupled (Wolpert et al. (1995), Huston and Jayaraman (2011)). The gustatory system is a perfect model to study this proposal, because animals necessarily respond to tastes with discriminative behavioral responses as they decide to swallow or expel the sensory stimulus in the mouth (Grill and Norgren (1978a), Katz and Sadacca (2011), ).
Previous studies of the rodent gustatory system have revealed that sensory-motor coupling is inherent in the temporal dynamics of gustatory cortical (GC) activity in response to taste stimuli. GC neurons respond to taste presentation with a sequence of firing-rate "epochs", two of which are taste-specific: neural firing in the first epoch carries information regarding the physio-chemical identity of the taste stimulus and in the second correlates with palatability, a variable intimately linked with the animal's decision to ingest or expel the taste (Katz et al. (2001), Fontanini and Katz (2006), Grossman et al. (2008), Piette et al. (2012), Sadacca et al. (2012)). Furthermore, ensemble analyses reveal that the transition between these two epochs happens suddenly and coherently across neurons (Jones et al. (2007), Sadacca et al. (2016)). The transition time, though highly variable in latency (between 0.5-1.5s post stimulus), is a strong predictor of the onset of the animal's consumption-related orofacial behavior (Sadacca et al. (2016)) -even when the timing of this behavior is manipulated by learning (Moran and Katz (2014)) or cueing ). GC neural ensembles appear to "hop" from one attractor state to another during taste processing, with the hop representing both the reaching of a swallow/expel decision and the emission of a motor signal to brainstem circuits that generate orofacial behavior (Miller and Katz (2010), Miller (2016)).
A direct prediction of this temporally dynamic model of gustatory processing, and most specifically of the suggestion that the transition between firing-rate epochs is of particular importance, is that well-timed perturbations of GC activity will affect the time course of a rat's taste-reactive ingestion-egestion processing. This prediction recently received indirect support when it was shown that optogenetic inhibition of the entire GC taste response ) modestly changes the probability of rejection behaviors in response to aversive tastes ("gapes", Grill and Norgren (1978a), ).
Such gross and broad perturbation of gustatory processing is an inadequate test of this very specific prediction, however: for one thing, multi-second inactivation likely has secondary effects that confound interpretation, particularly when dealing with an outcome variable (ability to gape) that is known to depend on an interconnected network of brain regions of which GC is a single node (Smith and St John (1999), Riley and King (2013), Samuelsen and Fontanini (2016)); in addition, it is impossible to disambiguate any epoch-or event-specific effects on consumption behavior using such lengthy manipulations of activity. A much more definitive test would involve using optogenetics to inhibit GC taste responses for short periods of time in specific processing epochs as awake rats processed and responded to a range of tastes.
Here we report the results of this precise experiment. We recorded the activity of GC ensembles while simultaneously inhibiting the firing of these neurons using optogenetics (specifically, with virus-delivered ArchT) for brief (0.5s) periods before, during or after the "hop" to the palatability-and decision-related attractor state. Our results provide strong support for the hypothesized importance of the transition time itself, and suggests that important pre-transition taste processing is performed in GC itself. These results provide a glimpse into the attractorlike dynamics underlying cortical taste processing -perturbations only have an impact before the system settles into the palatability and decision-related "stable" state. GC is likely just one participatory node in these attractor dynamics and it seems that behavioral control shifts to brainstem circuits once this stable state has been reached.

Subjects
Adult, female Long-Evans rats (n=5; 275-300g at time of virus injection; 300-350g at time of electrode implantation) served as subjects in our study (in our hands, female Long-Evans rats have proven more docile than males, but we have observed no sex differences in the basic cortical dynamics of taste responding). The rats were housed in individual cages in a temperature and humidity controlled environment under a 12:12h light:dark cycle. All rats were given ad libitum access to food and water before experiments started. Rats were weighed daily and observed to never drop below 80% of their pre-surgery weight. All experimental methods were in compliance with National Institutes of Health guidelines and were approved in advance by the Brandeis University Institutional Animal Care and Use Committee.
We also performed a set of control analyses of data taken from 10 adult, female Long-Evans rats, previously published in Sadacca et al. (2016) and .

Virus injections
We injected adeno-associated virus (AAV9) coding for ArchT and green fluorescent protein (AAV9-CAG-ArchT-GFP, 2.5 × 10 11 particles per mL) into the GC of rats. This AAV serotype has been shown to effectively spread to and infect all cell types (Aschauer et al. (2013)) in regions including GC (Maier et al. (2015); ). Rats were first anesthetized using a ketamine/xylazine mixture (1mL ketamine, 0.05 mL xylazine/kg body weight) delivered via an intra-peritoneal injection. Supplemental anesthetic injections were given as needed. The head was shaved, cleaned with an iodine solution and 70% ethanol, and positioned into the stereotax. We then excised the scalp and cleaned and leveled the top of the skull. Small craniotomies were drilled bilaterally over the location of GC (anteroposterior +1.4mm from bregma, mediolateral ±5mm from bregma; Paxinos and Watson (2007)), the meningeal tissues were gently excised, and virus was infused.
Prior to injecting the infusate, we suspended the virus particles in a solution of phosphatebuffered saline (PBS) and Oregon Green 488 (Invitrogen). We then lowered a glass micropipette (tip diameter: 10-20µm) filled with this solution into the centers of the craniotomies, and performed a sequence of 3 injections bilaterally into GC: at 4.9, 4.7 and 4.5mm ventral to dura; virus was injected in a series of discrete pulses (44 pulses/location, with 25nL per pulse, 7s between consecutive pulses = 1.1µL total volume injected per depth) controlled by a Nanoject III microinjector (Drummond Scientific). Following each unilateral set of injections, the micropipette remained in place for 5 min, after which it was smoothly removed over the course of 1 minute so that fluid would not spread back up the micro-pipette track. Craniotomies were sealed using a silicone sealant (Kwik-Sil, WPI).
Following injections, the scalp was sutured and the rat was given analgesic (meloxicam 0.04mg/kg), saline and antibiotic (Pro-Pen-G 150,000U/kg) injections. Similar antibiotic and analgesic injections were delivered 24 and 48 hours later. Rats were then allowed to recover for 4-6 weeks, in order to ensure adequate infection and subsequent expression of optical channels (ArchT) and GFP.
2.1.3 Opto-trode, intra-oral cannula and EMG electrode implantation 4-6 weeks after virus infusion surgery, rats were again anesthetized, and implanted with bilateral GC opto-trode bundles. Each bundle consisted of either 30 or 32 recording microwires (0.0015inch formvar-coated nichrome wire; AM Systems) and 1 optical fiber (0.22 numerical aperture, 200µm core, inserted through a 2.5mm multimode stainless-steel ferrule; Thorlabs). The microwire bundle was glued to a custom-made electrode-interface board (San Francisco Circuits) and connected to a 32 channel Omnetics connector. In the case of the 30 microwire bundles, the final two pins were connected to 2 electromyography (EMG) electrodes (PFAcoated stainless steel wire; AM Systems) implanted into the digastric muscle under the jaw. Finally, the microwires and optical fiber were connected to a custom-built 3D printed microdrive that allowed the entire assembly to be moved ventrally after implantation. The microwire tips were located 0.5mm ventral to the tip of the optical fiber -this maximized the likelihood that the electrodes recorded the activity of neurons that were illuminated by the laser. For more information on the implanted apparati and associated electronics, see Katz et al. (2001), Sadacca et al. (2016 and , as well as the Katz Lab webpage. The implantation surgery started similar to the virus injections: after anesthetizing the animal, we shaved and cleaned the scalp and situated the head in the stereotax. After excising the scalp and leveling the skull, we drilled 5 self-tapping screws into the skull for supporting and grounding the opto-trode bundles. The silicone seal was removed from the craniotomies, as were any tissues that had grown in since the prior surgery. We then slowly (over 5-10 minutes) lowered the opto-trode bundles to a depth of 4.3mm from the dura mater (0.2mm above the most dorsal location of virus injection). The ground wires were wound tightly around the skull screws and the bundles were cemented in place with dental acrylic. The optical fiber was looped so that the ferrule could be cemented away from the microdrive -this configuration reduced the stress on the microdrive when the animal was later plugged in to the experimental apparatus.
Once the opto-trode assembly was cemented in place, the rat was removed from the stereotax and implanted with a single (right-side) intra-oral cannula (IOC) for controlled delivery of tastants on the tongue. IOCs were made with thin polyethylene tubing and inserted in the space between the first maxillary molar and the lip, through the masseter muscle and inside the zygomatic arch, and out through the opening in the scalp (Phillips and Norgren (1970); Katz et al. (2001)) The IOC was topped with a plastic connector that could be attached to the taste delivery apparatus, and cemented in place with dental acrylic.
The EMG electrodes were channeled down the left side of the face (opposite from the IOC); after the overlying skin had been teased away from the belly of the digastric muscle, one end of each EMG electrode was tied to a suture needle, which was then inserted into the muscle, such that the electrode could be pulled into the desired position (for more details, see Loeb and Gans (1986); Travers and Norgren (1986); Dinardo and Travers (1994); ). The electrode wires were trimmed and held in place with vetbond tissue adhesive (3M) and the skin covering the anterior digastric was sutured back into place. Finally, a modified falcon tube was glued to the front of the headcap as a protective cap, and bacitracin ointment was applied all around the base of the headcap and over the wound under the jaw.
Following the surgeries, rats were given analgesic (Buprenophine 0.05mg/kg), saline, and antibiotic (Pro-Pen-G 150,000U/kg) injections. Similar antibiotic, saline and analgesic injections were delivered 24, 48 and 72 hours later, and bacitracin ointment was reapplied. The rats were handled every day and allowed to recover to 90% of their pre-surgery weight (at least 7 days after surgery) before being introduced to the experimental apparatus.

Habituation
Following recovery from the opto-trode implantation surgery, we habituated the rats to passive water deliveries for 3 days before any experiments started. In these daily habituation sessions, we attached the rats to the electrophysiology acquisition system, laser patch cables and taste delivery apparatus, and infused 100 pulses of distilled water (∼40µL per pulse; 15s inter-pulse interval) into the animal's oral cavity through the IOC. Starting with the second habituation day, we also placed rats on a mild water restriction schedule -20mL of water (not including the 4mL delivered during habituation sessions themselves) per day. This water restriction schedule was maintained for the duration of the study (∼7 days per animal).
Opto-trode bundles were driven deeper after each session using the microdrive built into the assembly; by the end of the habituation period, the distance traveled was 0.2mm, such that the tips of the electrodes lay within the region of GC infected with the virus.

Passive taste administration and laser stimulus delivery
We used 2 concentrations of palatable sucrose (30mM: Dilute Sucrose (Dil Suc), 300mM: Concentrated Sucrose (Conc Suc)) and of aversive quinine-HCl (0.1mM: Dilute Quinine-HCl (Dil Qui), 1mM: Concentrated Quinine-HCl (Conc Qui)) dissolved in distilled water as the stimuli in our experiments. Concentrated sucrose and quinine are rich in palatability-related valence and evoke strong orofacial responses; the dilute stimuli are of similar but less extreme palatability -a fact that aided in the analysis of palatability-related neural firing ; see also below). The taste delivery apparatus consisted of gently pressurized tubes containing taste solutions; the tubes converged upon a manifold of finer polyamide tubes that could be inserted into (to 0.5 mm past the end of) the IOC, thus eliminating any chance of mixing. The manifold could be locked securely into the dental acrylic cap. The tastes were then delivered under slight nitrogen pressure -this taste delivery protocol has been consistently shown to ensure reliable tongue coverage at short latencies (Katz et al. (2001);Sadacca et al. (2016); ).
We conducted 2 types of optogenetic perturbations in this study: 1) some sessions included trials with a "long" perturbation, with the laser turned on for the period of 0-2.5s post taste delivery; and 2) other sessions included trials with one of several "short" perturbations, for which the laser was turned on for 0.5s at either 0, 0.7 or 1.4s post taste delivery. We used a 532nm, DPSS laser (Laserglow Technologies), connected to the implanted ferrules using standard FC/PC patch cables (Thorlabs), for all optogenetic perturbations. The strength of the laser input was calibrated, prior to opto-trode implantation, to yield an illumination power of 40mW at the tip of the optical fiber. This output power perturbs all ArchT infected neurons in a 1mm 3 sphere below the tip of the fiber in vivo (Han et al. (2011);Yizhar et al. (2011)) -a sphere that encompasses about 33% of GC in the caudal-rostral axis (Kosar et al. (1986);Maier et al. (2015); ). These parameters have previously been shown to reduce the activity of ArchT+ cortical neurons with minimal latency and damage (Maier et al. (2015); Flores et al. (2018)).
We ran 1 experimental session per day per rat -only 1 type of optogenetic perturbation (either long or short) was done in each session. We interleaved the length of the optogenetic perturbation across sessions, so that a rat that experienced 2.5s perturbations in one session got 0.5s perturbations the following day. Taste and laser delivery were controlled through a Raspberry Pi computer.
Sessions with 2.5s perturbations consisted of 8 sets of trials (2 sets per taste -one with the lasers on and one with no laser). Each set included 15 trials, for a total of 120 trials per session. Similarly, sessions with 0.5s perturbations included 16 sets of trials (4 sets per taste -one with lasers on from 0.0-0.5s, one with lasers on from 0.7-1.2s, one with lasers on from 1.4-1.9s, and one with no lasers). To keep the total number of trials per session from ballooning (a basic concern in taste research is the awake animal's finite appetite), each set included only 8 trials (total, 128 trials per session). Again, we moved the opto-trode bundle 0.075mm ventrally (deeper into GC) prior to each session, to ensure that we obtained fresh units in every session.
Trials were delivered in pseudo-random order. Each involved delivery of ∼40µL of fluid through the IOC, for a total volume of 5mL per session.

Acquisition of electrophysiological data
We collected voltage samples from the implanted neural and EMG electrodes at 30kSamples/s using 32 channel analog-to-digital converter chips (RHD2132) from Intan Technologies. These chips are capable of recording voltage signals over a wide range of frequencies (0.1Hz-20kHz) and amplitudes (microvolts to millivolts), thereby enabling us to record neural and EMG signals through the same hardware system. The experimental chamber was ensconced in a Faraday cage that shielded recordings from external electrostatic and electromagnetic influences.

Histology and evaluation of GFP expression
In preparation for histology, rats were deeply anesthetized with an overdose of the ketamine/xylazine mixture, after which DC current (7µA for 7s) was passed through selected microwires, marking the area below the electrode tips. We perfused the rats through the heart with 0.9% saline followed by 10% formalin and harvested the brain. The brain tissue was incubated in a fixing mixture of 30% sucrose and 10% formalin for 7 days before GC was sectioned into 50µm coronal slices.
We rinsed the slices 3 times with 1X-PBS over 15 minutes and permeabilized them in a 0.3% Triton X-100+1% normal Donkey serum+1X-PBS blocking solution for 2 hours at room temperature. We replaced the blocking solution with primary antibody solution (1:500 anti-GFP-rabbit IgG; Life Technologies) for 12 hours at 4 • C. After incubation with the primary antibody, the slices were rinsed with 1X-PBS 3 times over 15 minutes followed by incubation with the secondary antibody incubation of (1:200 Alexa Flour 488 donkey anti-rabbit IgG (H+L); Life Technologies) for 12 hours at 4 • C. After a final set of rinses with 1X-PBS (3 times over 15 minutes), we mounted the slices on charged glass slides and cover-slipped them with Fluoromount Aqueous Mounting Medium. Slices were imaged with a Keyence fluorescence microscope to confirm successful virus infection and opto-trode location for each animal.
The spread of AAV in GC was evaluated via the expression of GFP, as done previously

Data analysis
Most statistical analyses were performed using Bayesian methods implemented in the PyMC3 probabilistic programming package (Salvatier et al. (2016)). The Bayesian framework entails the construction of flexible models of the data-generating process, and allows inference of the posterior distribution of the model parameters in light of the data actually observed (McElreath (2015)). We performed inference in these probabilistic models using the No-U-Turn-Sampler (NUTS; Hoffman and Gelman (2014)), a Markov Chain Monte Carlo (MCMC) algorithm that draws samples from the posterior efficiently. We confirmed the convergence of the sampling process by running multiple sampling chains per analysis and computing the Gelman-RubinR statistic (Gelman et al. (2011)) -R close to 1 (we allow values from 0.99 to 1.01) indicates that the sampling runs have converged and produced samples from the same posterior distribution. Finally, the results of the analyses are reported as 95% credible intervals for all inferred parameters. Credible intervals serve inherently as significance tests -for instance, a 95% credible interval that does not overlap 0 indicates statistical significance at the 5% level.

Single unit isolation
We designed a spike sorting toolbox in Python for electrophysiological data recorded using the Intan system. In brief, we followed a semi-supervised spike sorting strategy: the voltage data was filtered between 300-3000Hz and a Gaussian Mixture Model (GMM) identified potential clusters which were refined manually. For more details on our spike sorting methods and its efficacy in isolating single units, please consult Mukherjee et al. (2017). Our spike sorting code is freely available at blech_clust.

Impact of optogenetics on neural firing
We built a hierarchical Poisson generalized linear model (GLM) for the spiking of a single neuron to evaluate the impact of optogenetic perturbations on firing. Hierarchical GLMs provide precise estimates of condition-specific model parameters, especially when they are expected to vary around condition-agnostic means. In our case, the model parameters are the mean firing rates for every taste and optogenetic condition, that are in turn composed of taste-and optogeneticspecific effects ("random effects") and means across tastes and optogenetic conditions ("fixed effects"). Coupled with the Poisson distribution's suitability for count (here spikes) data, this model can accurately estimate the change in neural firing induced by optogenetic perturbations.
For each neuron n in our dataset, we aggregated its spikes on trial i of taste T in optogenetic condition O. There were 4 levels for T corresponding to the tastes in our dataset: Dil Suc, Conc Suc, Dil Qui and Conc Qui. The number of levels for O depended on the type of optogenetic perturbation being delivered in the session. In the 2.5s perturbation sessions, O had two levels, corresponding to the laser off (control) and on trials respectively. The 0.5s perturbation sessions, on the other hand, had 3 types of perturbations -starting at 0s, 0.7s or 1.4s after taste delivery. The 0.5s sessions, therefore, had 6 levels for O: a "laser off-laser on" pair for each of the 3 types of perturbations. Our model posits that the aggregate number of spikes S n,i,T,O is Poissondistributed with a mean (f iring n,T,O ) that depends on the taste (µ T ), optogenetic condition (µ O ) and an interaction between the taste and optogenetic condition (µ T,O ). As described above, owing to the hierarchical structure of the model, each of these effects is further composed of a fixed effect and a random effect. Using weakly informative Gaussian and Half-Cauchy priors for the mean and variance parameters respectively, our model formally says: Fixed effects: F 1 , F 2 , F 3 ∼ N (0, 10) Taste-and-optogenetics-specific means: As explained in the introduction to the data analysis section, we used MCMC to sample the posterior distribution of f iring n,T,O for every taste and optogenetics condition. We performed this analysis for every neuron in our dataset and finally calculated the impact of optogenetics on firing as the difference in f iring n,T,O between laser off (control) and their corresponding laser on trials. If the 95% confidence interval for these differences in f iring n,T,O for a neuron did not overlap 0, we concluded that the optogenetics significantly impacted the firing of this neuron.

Correlation of single neuron firing with palatability ranks
We analyzed the time course of palatability-related information in the activity of single neurons by regressing their firing rates on the palatability ranks of the tastes (Dil Suc: 3, Conc Suc:4, Dil Qui: 2, Conc Qui: 1; higher is more palatable). In order to estimate the firing rates of neurons, we aggregated the spikes of each neuron, on a trial-by-trial basis, in 250ms bins moved by 25ms steps. We divided the aggregate number of spikes by the width of the bins (250ms) to obtain the near-instantaneous firing rate of each neuron across time on individual trials.
These firing rates, of course, are likely to be very different between neurons. Furthermore, palatability-related firing may come in multiple forms -the correlation between firing rate and palatability ranks might be positive or negative (assuming that it isn't essentially 0) for different neurons, and be deemed similarly palatability-related. We therefore needed to perform a 2stage transform on neural firing before we could analyze the neurons together in our regression analysis. The first step was standardization -we transformed the firing rate of every neuron in each time bin by subtracting the trial-averaged firing rate of the neuron in that time bin and scaling by its standard deviation across trials (to get z-scores), ensuring that the firing rates of all neurons were on a comparable scale. Next, we multiplied the standardized firing rate of each neuron by the sign of the time-averaged Spearman correlation coefficient between its firing and the palatability ranks. This ensured that the sign of the relationship of neural firing with palatability was the same for all neurons in our dataset, but left the magnitude of that relationship unaffected.
Our statistical model treats the standardized firing rate f iring t,P,i of a neuron at time bin t on trial i of a taste with palatability rank P as Gaussian-distributed with a mean µ t,P that depends linearly on P . We defined the palatability index in time bin t, β P alatability,t , as the change in µ t,P induced by a unit change in P . β P alatability,t is, therefore, the slope of the line that explains µ t,P in terms of P , an estimate of the strength of the firing-palatability relationship. Using weakly informative Gaussian and Half-Cauchy priors for the mean and variance parameters respectively, our model formally says: Prior on palatability index: β P alatability,t ∼ N (0, 1)

Prior on observation noise: σ ∼ Half Cauchy(1)
Mean firing rate: µ t,P = β P alatability,t × P Firing rate: f iring t,P,i ∼ N (µ t,P , σ) We used MCMC to infer the posterior distribution of β P alatability,t across all neurons in our dataset (refer to the introduction of the data analysis section for more details on the use of Bayesian inference in our analysis). The firing rate transformations defined previously put the activity of all neurons on a comparable scale and allowed us to infer a single posterior distribution of β P alatability,t across all the neurons in our dataset. We repeated this regression for each time bin t from 0.25s before to 1.5s after taste delivery, obtaining posterior estimates of β P alatability,t specific to each time bin. Finally, we normalized β P alatability,t by subtracting its average baseline value (from 0.25 to 0s before tastes). We report the baseline-normalized β P alatability,t as the palatability index β P alatability .

Characterizing the time course of the palatability index
Similar to Sadacca et al. (2016), we modeled the time course of the posterior mean of the single neuron palatability firing index, β P alatability , with a logistic sigmoid. The difference between the lower and upper asymptotes of the S-shaped logistic function fits the total rise in β P alatability across time, while its slope describes the rate of this rise. As β P alatability was already normalized to its average pre-stimulus value, we set the lower asymptote of the logistic function to 0. With weakly informative Gaussian priors (restricted to positive values) on the upper asymptote (L), slope (k) and inflection time (t 0 , ms post taste delivery) of the logistic sigmoid, our model is as follows: We defined the peak of the palatability firing index, t peak , as the time (post taste delivery) when β P alatability reached 95% of its maximum value, L. We transformed the posterior distributions of L, k and t 0 to get t peak (inferred using MCMC) as follows:

Modeling and changepoint identification in ensemble firing data
As described in the Introduction (and Discussion), previous analyses reveal that GC population activity in response to a taste consists of a sequence of 3 coherent, suddenly-appearing ensemble states (Katz et al. (2001), Sadacca et al. (2016), ) in which firing rates code, in turn, for taste presence, taste identity, and taste palatability; the transition into this last state has particular relevance for the prediction of palatability-related behavior in single trials, and is the subject of this study. While identifying these sequences typically requires several passes through a dataset made up of many identical (i.e., unperturbed) trials, the good deal of work already done on the nature of these states (see also Jones et al. (2007) and Moran and Katz (2014)) renders it possible (for the purposes of the current study) to more concretely define this process as involving sudden coherent changes of ensemble firing into states, listed in their order of occurrence, in which the recorded ensemble has the following properties (also see Figure 6): 1. Detection state: The same distribution of population activity for all the tastes, indicating taste presence on the tongue.
2. Identity state: 2 distinct distributions of population activity, for the 2 taste identities in our experiments (Suc and Qui).
3. Palatability state: 4 distinct distributions of population activity, for the 4 taste palatabilities in our experiments (Dil Suc, Conc Suc, Dil Qui and Conc Qui).
With this characterization we were able to design a relatively simple changepoint model that allowed us to detect these coherent transitions in population activity in individual trials. We first prepared the data for the changepoint model by aggregating the spikes of each neuron in each trial into 10ms non-overlapping bins, indexing each neuron recorded in a session with a scalar i running from 0 to the number of neurons in the session N . We then converted the aggregate spiking data to a categorical format by marking each time bin by the index S of the neuron that spiked in that bin, with S = 0 corresponding to no spikes from any neuron. If more than one neuron spiked in a time bin -a highly uncommon occurrence, given the relatively low firing rates of GC neurons and the small (10 ms) bins being used -we randomly selected one of the spiking neurons for assignment to that bin (Jones et al. (2007); Sadacca et al. (2016)).
With the (processed) categorical spiking data in hand, we now designed a changepoint model to describe the ensemble firing in each of the 3 states (listed above) as categorical distributions with N +1 emissions, with 1, 2 and 4 such distributions corresponding to the detection, identity and palatability states respectively. We analyzed 1.5s of ensemble activity post taste delivery from each of the 4 optogenetic conditions in the 0.5s perturbation sessions. For the control (laser off) trials, this corresponded to 0-1.5s of firing after taste delivery. On the perturbed trials, we ignored the 0.5s of activity when the lasers were on -for example, we analyzed 0.5-2.0s of firing post tastes when the lasers were on from 0-0.5s. In the resultant 1.5s of activity, we assumed that the change from detection to the identity state, C I , happens anywhere in the interval [0.2s, 0.6s] (except the 0-0.5s perturbation trials, where the identity state can start earlier at 0.1s). The second changepoint, C P , from identity to palatability firing was assumed to occur anywhere in the interval [C I +0.2s, 1.3s] (except the 0.7-1.2s perturbation trials, where the palatability state can start earlier at C I + 0.1s). This is equivalent to placing uniform priors over the intervals that define C I and C P , and corresponds to the timing of sudden, coherent firing rate transitions in GC ensembles (Jones et al. (2007), Sadacca et al. (2016)).
C I and C P are, therefore, latent variables of the changepoint model that control the probabilities of the emissions actually observed. The Expectation-Maximization (EM) algorithm is a widely used approach to perform inference in such models with latent variables -for stability and speed issues, however, we used a "hard-assignment" version of EM to fit the changepoint model (Bishop (2016)). Starting with a randomly chosen set of initial emission probabilities α D , α I and α P for the categorical emissions that define the detection, identity and palatability states respectively, the EM algorithm for our changepoint model repeatedly cycled between 2 steps: 1. "Hard" E-step: Pick the combination of the latent variables, C I and C P , that has maximum posterior probability given the observed categorical spikes S and the ensemble firing probabilities α D , α I and α P . We directly pick the mode of the joint posterior distribution of C I and C P in our hard-assignment version of the E-step instead of taking their expectations/means.
2. M-step: Set the categorical firing probabilities for each state to values that maximize the likelihood of the data given the (C I , C P ) pair picked in the E-step. This is proportional to the number of emissions of each neuron in that state. For example, with S t as the emission observed at time t, the likelihood-maximizing emission probabilities of neuron n can be calculated as: In identity state: In palatability state: α P,n = 1.5s where 1 is the unit function that is 1 when S t = n and 0 otherwise.
In order to deal with the possibility that EM can get stuck at sub-optimal local maxima of log likelihood, we ran the algorithm from 100 different random initializations of the α parameters. We monitored the log likelihood of the data given the model parameters and ran the algorithm to a convergence threshold of 10 −8 (or a maximum of 300 iterations). Finally, we picked the run with the maximum log likelihood at convergence and reported the changepoints (and their posterior probabilities given S and α) found on this run.

Measuring aversive orofacial behaviors (gapes)
Bitter (e.g., Quinine) tastes cause rats to produce an orofacial behavior known as "gaping", the purpose of which is to maneuver the offending substances to the front of the mouth for possible egestion. As such, gapes index the fact that the neural processing of the bitter taste has in a certain sense reached completion -the rat has "decided" that it does not want to ingest the taste. The occurrence of gapes can be measured in a number of ways, the most common of which is via human coding of video recordings -in the best of circumstances, gapes are readily visible as large yawn-like movements.
Of course, the best of circumstances often fail to occur in rats free to move and rear. This fact, and the difficulty of getting precise measures of gape onset time from a visual record, renders video coding of gapes suboptimal for our purposes. Much more objective and less noiseridden is evaluation of jaw electromyography (EMG), in which individual gapes are recognizable as particularly large-amplitude and large-duration electrical bursts ( Figure 4A1-A2). We have previously built a quadratic classifier to detect these bursts in ongoing anterior digastric EMG signals, achieving 75% accuracy ).
Even this approach has somewhat troubling limitations, however, as its limited accuracy indicates. This limited accuracy stems from the facts that: 1) not all high-amplitude jaw movements are gapes; and 2) gapes vary widely in amplitude, and in fact some are small enough to appear similar in size to many other mouth movements (see Figure 4A1-A2). In practice, both types of variability leave the classifier subject to false positives that must be somehow recognized and removed -the former most notably at the beginning of trials (when the taste hits the tongue, causing 1-2 relatively large-amplitude licks), and the latter in responses to Sucrose (among other situations).
One solution to these problems involves making simultaneous recordings from multiple jaw muscles, but pilot experiments left us concerned that such drastic infiltration of the jaw can compromise normal movement, which would make interpreting our results difficult. Instead, we decided to take advantage of another, more robust feature of gaping: the fact that gapes occur in 4-6 Hz "bouts" of anterior digastric activity (Travers and Norgren (1986), ). While identifying gaping bouts as time periods during which this rhythmicity dominates the EMG signal is also imperfect -it is probabilistic and involves smoothing across time -it solves most of the problems described above.
We instantiated just such an procedure here, applying a Bayesian spectrum analysis that estimates the posterior probability that a 4-6Hz rhythm underlies a short time series of EMG activity (see below for technical details). By this analysis, the probability of gaping to any taste is modestly elevated at trial onset (because of the initial large-amplitude licks), but it quickly drops to effectively zero for Sucrose, contributing nothing to the overall calculation of when gaping begins. On Quinine trials, in contrast, the probability accords well with gape bouts ( Figure 4B1-B2), rising precipitously and reliably just prior to the first gape (detected in a subset of data with both video recordings and the quadratic classifier, Figure 4D).
Finally, as both the probability and onset of gaping behavior differs between Dil and Conc Qui (Grill and Norgren (1978a), Travers and Norgren (1986)), we defined the overall average gaping latency as the difference between the two individual time distributions of gaping probability, which can be statistically assessed as the Kullback-Leibler (KL) divergence (again, see technical details below). Not only does this procedure calculate the onset of aversive orofacial behavior in terms of when the two (differently) aversive tastes in our battery can be distinguished on EMG, it provides a natural threshold for this comparison. By pitting the two Qui concentrations against each other, this method gets rid of most of the nonspecific gape-like EMG activity (mentioned above) which is of similar magnitude on both Dil and Conc Qui trials and does not contribute to the gape onset calculation. Contrary to previous methods, which (usually) removed trials where gapes could not be reliably detected from further analysis, this algorithm combines EMG data from all the trials available and avoids statistical comparisons between conditions with very different sample sizes. Thus, at the cost of being unable to precisely detect the specific timing of later gapes in a bout, this procedure provides a robust estimate of the average timing of the first gape.
Bayesian spectrum analysis (BSA) of EMG recordings As detailed previously, we recorded voltage signals from 2 unipolar EMG electrodes implanted in the anterior digastric muscle at 30kSamples/s. We used the difference in the voltage recorded by the 2 electrodes as the EMG signal -this procedure helps to cancel out any large artifacts produced by the animal's movements and is equivalent to using a differential amplifier (as done in ). We down-sampled the EMG signal to 1000Hz by averaging the voltage values in sets of 30, and highpass filtered the down-sampled signal above 300Hz (Travers and Norgren (1986); ) using a 2 nd order Butterworth filter. The absolute value/magnitude of the filtered EMG signal was then lowpass filtered (again using a Butterworth filter of order 2) below 15Hz, effectively capturing the envelope of variation of the EMG signal. This cutoff of 15Hz is sufficient for identifying orofacial behaviors, all of which occur at frequencies smaller than 10Hz (Grill and Norgren (1978a); ).
We subjected the envelope of the EMG signal to Bayesian spectrum analysis (BSA). BSA involves the construction of a probabilistic model of the generation of periodic signals from the superposition of sinusoids of different frequencies. We divided the signal on each trial into bins of width 300ms, with a step size of 1ms. We assumed that the EMG signal in each bin is produced by a sinusoid of a single frequency (plus noise) -in a probabilistic setting, this assumption implies the same model as a discrete-time Fourier transform. Contrary to the Fourier transform, however, BSA infers the posterior distribution of frequencies given the data. BSA has been shown to provide posterior estimates of frequencies that are an order of magnitude more precise than the Fourier transform (Bretthorst (2013); Granqvist et al. (2011)). We used the BaSAR R package for BSA (Granqvist et al. (2012)) and calculated the posterior probabilities of frequencies from 1Hz to 10Hz in 20 steps for each 300ms wide bin of data.
Identifying the mean onset of aversive orofacial behavior Rats respond to intra-oral deliveries of Qui in the concentration range used in our experiments (10 −4 to 10 −3 M) with an initial set of non-specific investigative licks that are followed by large, jaw-opening mouth movements called gapes (Grill and Norgren (1978a), Figure 4A1-A2). Gapes primarily involve activity of the anterior digastric muscle at 4-6Hz (Grill and Norgren (1978a), ) -we, therefore, used the probability of movements at 4-6Hz in the digastric EMG signal (from BSA, see previous section) as the probability of gaping (Pr Gape ). This spectral measure of Pr Gape has a strong correspondence with a previously-defined (and above-discussed quadratic classifier that tags individual mouth movements as gapes ). On individual Qui trials ( Figure 4B1-B2), Pr Gape from BSA is high (close to 1.0) when the quadratic classifier tags mouth movements as gapes. In addition, the average probability of gaping (Pr Gape ) from BSA ( Figure 4C1-C2) is very similar to an across-trial, peri-stimulus average of the gapes picked by the quadratic classifier. In contrast to the quadratic classifier, however, the BSA measure of Pr Gape is based entirely on the spectral content of the EMG signal. It, therefore, does not require the construction of a sufficiently complex classifier function (with a large enough set of experimenter-tagged examples to train the classifier) to pick out individual gapes. This also ensures that BSA considers bouts of movements together while calculating Pr Gape , making it robust against isolated large amplitude movements early in the animal's orofacial response. These initial movements were often found to be large licks on video and limited the accuracy of the quadratic classifier in  to 75%.
The probability of the transition from the rats' initial investigative licks to gapes depends on the concentration of Qui delivered: 10 −3 M (Conc Qui) elicits gapes on more than twice the number of trials as 10 −4 M (Dil Qui) (Grill and Norgren (1978a), ). Comparison of Pr Gape on Dil and Conc Qui trials, thus, provides a natural way to calculate the mean onset of gaping across all the Qui trials in an experimental condition. We expect the distribution of Pr Gape on Dil Qui trials to be similar to that on Conc Qui trials in the investigative licking phase. Once gaping starts, however, we expect a large difference in the distributions of Pr Gape on Dil and Conc Qui trials. Pr Gape on Dil Qui trials, therefore, acts like a baseline for Pr Gape on Conc Qui trials: we conclude that gapes have started only when Pr Gape of Conc Qui begins to differ significantly from this baseline.
We used Beta distributions to describe Pr Gape on Dil and Conc Qui trials. The Beta distribution is commonly used to model the probability parameter of a Bernoulli (1/0) process 1 . Gaping being a Bernoulli process, the Beta distribution is an appropriate choice for modeling Pr Gape . We defined one such Beta distribution in each time bin for Dil and Conc Qui separately, parametrized by the number of trials where the animal was gaping (Pr Gape > 0.5) or not (Pr Gape < 0.5). The Kullback-Leibler divergence of these Beta distributions (D KL (Conc Qui||Dil Qui)) 2 provides a natural way to quantify the difference between Pr Gape on Dil and Conc Qui trials and shows a sharp jump ∼1s post taste delivery ( Figure 4E), consistent with the timing of the transition from investigative licks to gapes (Grill and Norgren (1978a), Travers and Norgren (1986), ). Finally, we calculated the cumulative sum of D KL (Conc Qui||Dil Qui) across time: the jump corresponding to the mean onset of gaping is expressed as a change in slope of the cumulative sum. We fit two straight lines to the cumulative sum to capture this change in slope: the intersection of the two lines defines the mean timing of the onset of gaping ( Figure 4F).
1 The Beta distribution for the parameter p of a Bernoulli process is expressed in terms of its concentration parameters, α and β. α = observed number of 1s and β = observed number of 0s.
2 The KL divergence between two Beta distributions with concentration parameters (α 1 , β 1 ) and (α 2 , β 2 ) can be written as: where Γ and ψ are the gamma and digamma functions respectively. Figure 1A depicts the preparation used for our experients -IOCs for taste delivery, bilateral GC opto-trodes for recording of neural ensemble activity and delivery of laser light, and EMG electrodes in the anterior digastric (jaw) muscle for simultaneous assaying of consumption-related mouth movements. Four weeks prior to the surgery in which we installed these assemblies, we had injected AAV carrying the optogenetic silencer ArchT (along with green fluorescent protein -GFP) into GC. The latter allowed us to confirm (post-mortem) the expression of ArchT in GC neurons by immunohistochemical verification of the GFP tag ( Figure 1B). The rats received intra-oral deliveries of 30mM sucrose (Dil Suc), 300mM sucrose (Conc Suc), 0.1mM Quinine-HCl (Dil Qui) and 1mM Quinine-HCl (Conc Qui). On 3 4 of the deliveries (i.e., trials), we inhibited GC neurons for 0.5s, beginning either at 0s, 0.7s or 1.4s post taste delivery ( Figure 1C). These three perturbation windows tile the period containing the temporal epochs that characterize GC taste responses (Katz et al. (2001) 2016)): more specifically, the earliest (0-0.5s) and latest (1.4-1.9s) inhibitions affect GC neurons before and after the range of likely transition times into palatability-related firing, which typically occur just before or during the middle (0.7-1.2s) perturbations. In a separate set of experimental sessions (performed using the same rats), we inhibited GC across the entire duration of the taste responses (0-2.5s post stimulus) ( Figure 1C) as a control comparison for the brief 0.5s perturbations.

Experimental paradigm and data overview
We recorded the activity of 244 GC single units across 10 experimental sessions (24.4 ± 13 units/session) of 0.5s inhibition. We recorded the activity of an additional 73 GC single units in 5 experimental sessions (14.6 ± 4.7 units/session) of 2.5s inhibition. The two types of experimental sessions were counterbalanced, such that 3 of the 5 rats received 2.5s inhibition sessions first, and 2 received 0.5s inhibition sessions first. No differences with order were noted. The AAV-ArchT construct used in this study has been shown to infect neurons of multiple types (e.g., pyramidal neurons and interneurons) in an unbiased manner (Aschauer et al. (2013)). Our optogenetic inhibition protocol, therefore, can be thought of as a general perturbation of the dynamics of GC neurons in response to tastes -a general perturbation that would be expected (perhaps paradoxically) to enhance the firing of some neurons through networklevel effects (like disinhibition, via suppression of the firing of inhibitory neurons, Allen et al. (2015)). This expectation was borne out in the data: the firing of most of the recorded GC units (146/244, 60%, example unit in Figure 2A1-A4) was significantly suppressed when the laser was switched on for 0.5s, but the firing of an additional 20% (49/244) was significantly enhanced.
The same pattern of results was observed when the duration of optogenetic inactivation was increased to 2.5s: the firing of 82% of GC neurons (60/73, example unit in Figure 2B1-B2) was inhibited, and the activity of 15% (11/73) was enhanced. The fact that 2.5s of laser stimulation appeared to inhibit a larger percentage of neurons is likely an artifact of analysis methods: suppression of the low firing-rates (3-10Hz) that dominate GC taste responses (Katz et al. (2001), Jones et al. (2007)) can be difficult to detect, particularly in short windows; consistent with this, we observed that the highest likelihood of detecting suppression in 0.5s perturbation sessions occurred when that perturbation was delivered in the middle of taste processing (0.7-1.2s, Figure 2C) -at the time of peak firing rate modulations. With 2.5s of inactivation, which covered the entirety of GC taste responses, we naturally had the power to detect suppression in a larger fraction of neurons ( Figure 2D).
Although this specific optogenetic protocol cannot be used to answer cell-type/microcircuitspecific questions, its network-wide effects are ideal for testing the dynamical systems nature of taste processing in GC (the purpose of the current work): GC taste responses evolve through a sequence of temporal epochs (Katz et al. (2001)) which have the hallmarks of emergent, quasistable states of a system that can be speculatively described, at a high level, as an attractor network (Jones et al. (2007), Miller and Katz (2010), Sadacca et al. (2016)). Our optogenetic protocol brings about a massive perturbation of the network activity characterizing these stable states; by mapping the state dependence of the effects of these perturbations, we are able to directly test the proposed existence of these attractor states.

Early perturbations delay single neuron palatability responses while late perturbations do not
We first assessed the impact of optogenetic perturbation on the palatability-related content of GC taste responses that had been smoothed (using 250ms-wide windows moved in 25ms steps) and standardized to be on a uniform scale (see Materials and Methods for details). The set of responses (1 per taste) were regressed against the palatability ranks of the taste stimuli (Conc Suc:1, Dil Suc:2, Dil Qui:3, Conc Qui:4) to obtain a palatability index, β P alatability . This Bayesian regression analysis gives access to the entire posterior distribution of β P alatability at every time point; we can, therefore, conclude that β P alatability is significantly different from 0 if the 95% credible interval of its posterior distribution does not overlap 0 (marked by dots in Figure 3A). We used logistic sigmoid functions to better characterize the time evolution of the posterior mean of β P alatability (shown with dashed lines in Figure 3A). We defined the size and latency (time to attain 95% of maximum size) of the upper asymptote of the logistic fit as the magnitude and latency of the peak of β P alatability respectively.
As expected, the 2.5s perturbation had devastating impact on palatability responses of neurons in the affected GC network ( Figure 3A). In line with previous studies (Sadacca et al. (2016)), β P alatability climbed to an asymptote ∼0.8s after taste delivery when the lasers remained off. However, on trials where the lasers were switched on at the time of taste delivery and left on for 2.5s, β P alatability never rose significantly from 0. Note that the latency to peak palatability firing seems comparable in the two conditions (blue bars in Figure 3B), but that the magnitude of the peak is effectively 0 when GC neurons are being perturbed (red bars in Figure 3B). The long perturbation, therefore, appears to impact the entirety of the GC taste response.
The impact of brief (0.5s) perturbations on the palatability content of GC taste responses was smaller in magnitude, but could be quite dramatic with regard to peak timing, depending on when the perturbation occurred ( Figure 3C). In these sessions, just as in the 2.5s perturbation sessions, β P alatability peaked ∼0.8s after taste delivery when the lasers were left off. Furthermore, neither the timing nor the magnitude of this peak were significantly affected by perturbation of GC neurons in the later part of the taste response (1.4-1.9s, after palatability-related firing had mostly emerged). In contrast, if activity was perturbed for the first 0.5s of the GC taste response, the palatability content of this response did not reach asymptote until ∼1.3s, a lag of almost 0.5s compared to the control condition (laser-off trials). Surprisingly, despite delaying the peak of β P alatability , the early perturbation did not affect its later emergence -if anything, the magnitude of the peak was larger in this condition (red bars in Figure 3C). The 0-0.5s perturbation, thus, appears to produce a transient shift out of the attractor dynamics of the GC taste response followed by relaxation back into the stable state at the end of the perturbation. The variability in the relaxation process itself (like overshooting the stable point, depending on the speed of relaxation) could easily explain the apparent increase in the magnitude of the peak palatability index in this condition. Finally, 0.5s perturbations delivered in the middle of the taste response (0.7-1.2s) had a powerful impact on GC palatability-related firing: the magnitude of the peak of β P alatability was significantly lower in this condition (red bars in Figure 3C); the latency of this peak, meanwhile, was (like that produced by earlier perturbations) about 0.5s later than no-laser trials. The former effect was unsurprising, as this particular perturbation overlaps the heart of palatability related activity in GC neurons (Katz et al. (2001), Sadacca et al. (2016). The impact of this GC perturbation was short-lived, however, and did not prevent the palatability content of GC responses to rise to a late asymptote once the lasers were switched off.

GC perturbation delays the onset of aversive orofacial behavior
We monitored our rats' mouth movements via electromyography (EMG). Specifically, we implanted EMG electrodes in the anterior digastric muscle; as a jaw moving muscle, the anterior digastric plays a major role in driving "gapes", the rhythmic orofacial behavior that serves to move aversive tastants to the front of the mouth in preparation for expelling. Far less accessible tongue muscles underlie mouth movements that support behaviors (such as "lateral tongue protrusions") that help the rat prepare to ingest appetitive tastants (Grill and Norgren (1978a), Travers and Norgren (1986), ). For that reason, we focus solely on gapes in this study.
Individual mouth movements can be recognized as bursts of anterior digastric EMG activity ( Figure 4A1-A2). However, the variability in the amplitudes and durations of these EMG bursts limits their usefulness in separating gapes from other large mouth movements. We, therefore, made use of a more robust feature of gaping -the fact that gapes occur in 4-6Hz bouts of several large mouth movements (Travers and Norgren (1986), ). We analyzed the spectral content of the envelope of the EMG signal using Bayesian spectrum analysis (BSA; see Materials and Methods for a detailed discussion) and measured the probability of gaping as the total posterior probability of 4-6Hz movements.
On individual Qui trials ( Figure 4B1-B2), this estimate of the probability of gaping has strong correspondence with the onset latency of gaping bouts identified by a classifier trained on individual bursts of EMG activity ); the trial-averaged probability of gaping calculated by BSA and the classifier are also similar, for both trial types in which gaping occurred (Dil and Conc Qui trials, Figure 4C1-C2). Finally, the fact that the probability of gaping jumps precipitously right before the first gape (as identified both on video and by burstbased classification, Figure 4D) confirms this algorithm's reliability in identifying periods of gaping in the EMG signal (see Materials and Methods for more details).
With this information in hand, we were able to investigate the effects that perturbations of GC activity have on the animals' rejection of aversive Qui, first showing that gaping on average begins ∼0.9 sec after Qui delivery, when analysis is restricted to the 25% of trials in which the laser was off (i.e., trials in which GC neurons were not perturbed, (Figure 5)A). This latency is consonant with a great deal of research using video analysis (Grill and Norgren (1978a)) or classic burst-oriented analysis of EMG (Travers and Norgren (1986)). Furthermore, this estimate matches that observed in control rats (published in Sadacca et al. (2016) and ) that received neither laser nor ArchT expression -thus we are able to conclude that, at least with regard to the driving of gaping, 0.5 sec of optogenetic inactivation does not have an impact on GC that alters later trials.
While previous work has shown that this onset of gaping closely follows the appearance of palatability-related firing in GC (which happens suddenly and coherently across neurons in single trials, see below and Sadacca et al. (2016)), GC does not seem involved in controlling the mechanics of individual gapes when they occur (Grill and Norgren (1978b), ). We therefore predicted that GC perturbations occurring well after the time of palatability appearance would have minimal impact on gaping behavior. In fact, our data show that rats gaped normally, with gape bouts beginning at approximately the same time as in control (no laser) trials, if perturbations arrived late in the trial (1.4-1.9s, Figure 5B). Further analysis (not shown) confirmed that this late perturbation failed to end gaping bouts that had already begun. GC, therefore, appears not to be required for the maintenance of gaping.
In contrast, GC plays a clear role in the initiation of gaping. GC perturbations timed to occur before (0-0.5s) the emergence of palatability information in GC activity, for instance, delayed gaping onset by approximately 0.25s on average ( Figure 5B). This delay cannot be explained in terms of removal of early gaping (which seldom occurred as early as 0.5s after taste delivery) -an analysis of control (no laser) trials showed that removing latencies of less than 0.5s did not change the mean onset time of gaping significantly. The much more likely explanation is that GC inhibition (which is inevitably partial, see Discussion) perturbs the ongoing process that leads to the release of a "decision to gape" signal visible in GC (Sadacca et al. (2016)); these data therefore confirm the involvement of GC in this process.
Similarly, GC perturbations timed to occur squarely around the time of the neural state change associated with behavior (0.7-1.2s; see Sadacca et al. (2016)) delayed the onset of gaping until just less than 1.2s after taste administration -approximately 0.25s after gaping on control trials and in control sessions. Not only was this impact of brief optogenetic perturbation significant, it was every bit as large as that observed with whole-trial (i.e., 2.5s) perturbations, which delayed the appearance of gaping by ∼0.2s without significantly altering the overall fraction of trials on which gapes occurred (data not shown). These long perturbations are not discussed further, because they had the additional unintended consequence of impacting gaping behavior on control trials (see Figure 5A and Discussion). That is, a range of brief disruptions of GC activity occurring in the first 1.2s of quinine processing had strong and similar impacts on the triggering of aversive orofacial behavior.

GC perturbation affected orofacial behavior only if it was delivered before the onset of palatability-related ensemble activity
We have previously demonstrated that changes in firing rates of GC neurons -changes that appear gradual in trial-averaged analysis -occur as sudden transitions between (primarily two) ensemble firing rate "states" (Jones et al. (2007)). Furthermore, not only are the temporal dynamics of GC taste responses better described in terms of transitions between these two states; in addition, the onset time of the latter state, which is laden with information about stimulus palatability and appears, on an average, around the time when single neuron palatabilityrelated firing peaks (∼1s post stimulus, Figure 3), is highly predictive of the latency of gaping on single trials (Sadacca et al. (2016)). Trial averaging smears the firing data into a gradual "ramp" because these transitions happen at wildly different latencies on different trials.
We timed our 0.7-1.2s perturbations to overlap with the transition into this palatabilityrelated ensemble activity state; due to the afore-described variable timing of this transition, however, there are likely to be a subset of trials in which the ensemble state emerged before the perturbation started at 0.7s. This fact affords us an opportunity: we predicted that the 0.7-1.2s perturbation would impact gaping latency differently depending on whether the transition into the late ensemble activity state had already occurred in that trial (and, by extension, that the trial-averaged results in Figure 5B occlude our ability to see multiple effects).
While we have previously used Hidden Markov Models (HMMs) to detect ensemble activity patterns in GC population responses to tastes (Jones et al. (2007), Katz (2014), Sadacca et al. (2016)), this analysis is not amenable to the data in our study -the optogenetic perturbations mean that HMMs do not have access to uninterrupted sequences of GC activity and cannot converge to stable fits. Instead, we took advantage of the insights gained from our previous publications (Katz et al. (2001), Katz (2006), Jones et al. (2007), Grossman et al. (2008)) to help us build a change-point model of GC population activity; specifically, the model consisted of 2 activity change-points, the latter of which led into palatability-related firing. This model constrained the general HMM framework in a way that allowed us to estimate whether the transition into the behavior-related state had occurred prior to the onset of optogenetic perturbation (i.e., 0.7s, see Figure 6 and Materials and Methods for details).
The distributions of transition times (identified by the change-point model) into the palatability-related ensemble state when GC neurons were perturbed from 0.7s to 1.2s post stimulus are shown in Figure 7A for all Qui trials. As firing rates were suppressed during the perturbation, we did not attempt to identify change-points when the lasers were on. The model found that the palatability state emerged before the lasers were switched on at 0.7s on 55% of the trials (92/168); on the remaining 45% of trials (76/168), the palatability change-point could be identified only after the laser was switched off at 1.2s. We confirmed that the early observed change points were in fact transitions into the palatability-related ensemble state: regression analyses revealed significant palatability information present before 0.7s in trial-averaged single neuron firing ( Figure 7B) for trials in which the ensemble state transition occurred prior to laser onset time; this information was notably lacking in the remaining trials. Finally, we divided these trials into two groups on the basis of the timing of this neural population transition and found, in line with our expectations, that identical 0.7-1.2s perturbations had distinctly different effects on the onset of gaping for the two groups of trials ( Figure 7C). Perturbations that arrived before the ensemble transition delayed the onset of gaping by more than ∼0.5s -that is, the average onset of gaping across this group of trials was > 0.2s after the end of GC inhibition. This does not simply represent a truncation of the control distribution of gaping latencies: we looked at gapes in control trials (both within brief inactivation sessions and in sessions with no laser); even when we restricted ourselves to studying only the trials in which gaping latency was later than 1.2s (the laser off time), the average gaping latency (∼1.3s, i.e, about 0.4s later than no-laser trials) was still significantly less then that observed in the subset of trials in which the ensemble transition failed to precede 0.7-1.2 GC inhibition.
In fact, only about 20% of no-laser trials lacked gaping-related EMG activity till as late as the average latency in this particular group of laser-on trials.
Gaping in trials in which the ensemble transition to the high-palatability state preceded the onset of GC inhibition at 0.7 occurred significantly earlier ( Figure 7C). Contrary to our expectation, however, gaping was somewhat delayed compared to the no-laser condition even in these trials. As the ensemble state transition necessarily happens by 0.7s on these trials (i.e., far earlier than the average transition time on control trials), we expected that the onset of gaping would be similarly expedited. This was not the result that we obtained. We considered several possible explanations for this result (see Discussion), the most reasonable of which seemed the possibility that ensemble transitions are likely non-instantaneous (and/or that their onsets are only noisily estimated by our change-point algorithm); if so, it is possible that transitions identified as occurring at (or just before) 0.7s might actually not be completed until after GC inhibition has begun, and that gaping on those trials is substantially delayed.
We tested this hypothesis, and found that indeed the delay in the onset of gaping can be entirely attributed to the trials where the ensemble state transition is calculated to occur between 0.65 and 0.7s. That is, the onset of gaping in trials in which the ensemble transition happened at 0.65s or earlier occurred more than 300ms earlier than in control trials. This result clearly shows that the onset of palatability-related population activity in GC marks a discrete shift in taste processing; if the laser perturbation begins even 50ms after this state transition, it fails to impact the timing of aversive orofacial responses.
Thus, our results, taken together, demonstrate not only that the result of brief optogenetic inhibition of GC depends on when within a single trial that inhibition occurs. In addition, they show that: 1) the ensemble transition in taste-related firing that predicts behavior is in fact the specific time at which the decision to gape is released; and 2) given the trial-to-trial variability in the issuance of this decision, the result of reliably-timed brief optogenetic inhibition of GC will vary from trial to trial -the impact of neural perturbation depends on precisely what state the brain has achieved prior to that perturbation.

Discussion
Perception and action are inextricably linked in cortical taste responses. Neurons in gustatory cortex (GC), the primary sensory cortical area for taste, exhibit temporally dynamic responses that gradually (across 2.5s post stimulus) shift from reflecting stimulus identity to predicting a rat's consumption decision in trial-averaged analyses (Katz et al. (2001), Fontanini and Katz (2006)). Upon being subjected to ensemble analysis, these gradual changes in firing rate epochs are revealed to be discrete and coherent transitions between population activity "states" (Jones et al. (2007)) -transitions that are so variably timed that stimulus-aligned averages effectively blur them out. Despite their highly variable latencies, these ensemble firing states reliably precede the onset of ingestion-egestion mouth movements by ∼0. 2-0.3s (Sadacca et al. (2016), ) -predicting not only the nature but the latency of these movements in single trials.
Here we show that the variably timed activity states of GC neural ensembles are not merely "efferent copy" reflections of processes occurring elsewhere, but instead indicate processing that is intrinsic to GC. We find that brief (0.5s) optogenetic perturbations of GC neurons impact the timing of the animal's decision to ingest or expel the taste in the mouth -but only if the perturbations begin before the neural ensemble has shifted to palatability-related firing. Thus, a unique moment in time (the shift of population activity to reflect stimulus palatability) that is enormously variable from trial-to-trial reflects a tipping point in taste processing; cortical disruptions have no impact beyond this tipping point, as the control of the ongoing movements themselves shifts elsewhere (presumably to brainstem pattern-generators that produce ingestion-egestion mouth movements soon after).
A massively interconnected network of brain regions underlies or reflects taste processing -apart from GC, this network includes the central and basolateral nuclei of the amygdala (CeA and BLA, Nishijo et al. (1998), Grossman et al. (2008), Fontanini et al. (2009), Sadacca et al. (2012), hippocampus (Stone et al. (2005), Ho et al. (2011), lateral hypothalamus (LH, Yamamoto et al. (1989), Li et al. (2013)), the nucleus of the solitary tract (NTS, Di Lorenzo and Lemon (2000)) and bed nucleus of the stria terminalis (BNST, Norgren (1976), Li and Cho (2006)). Several of these brain regions send direct descending feedback to the brainstem, influencing both its activity (Di Lorenzo (2000), , Cho et al. (2003, Li et al. (2005), Baez-Santiago et al. (2016)) and generation of orofacial movements (Zhang andSasamoto (1990), Shammah-Lagnado et al. (1992), Travers et al. (1997), Berridge and Valenstein (1991)). Given this widely distributed network of processing nodes, it is to be expected that perturbation (or disruption over long periods of time) of one (or a few) of the participatory nodes will initiate homeostatic mechanisms that guard against the degradation of behavior; thus, it is unsurprising that ablation (King et al. (2015)) or disruption of GC ) has, at best, modest effects on orofacial behavior -in fact, the basic gaping response to quinine is produced even in decerebrate rats (Grill and Norgren (1978b)). Nonetheless, we find that brief perturbations of GC do significantly alter these behaviors, proving that far more than the minimal circuit is involved in triggering them in situ.
Longer disruptions of GC activity appear to have spillover effects that can confound the interpretation of their behavioral impact -our 2.5s long optogenetic perturbations delayed the onset of gaping even on control (no laser) trials. Such spillover effects may reflect cellular or network-level processes, but they cannot be attributed to cell death caused by the perturbation: in our case, this optogenetic protocol has been shown to have no observable impact on cell integrity in GC, even for perturbations much longer than 2.5s (Maier et al. (2015), Flores et al. (2018); furthermore, the same rats in later sessions produce normally-timed orofacial responses on the control trials, despite the fact that gaping latencies were as altered by 0.5s GC perturbations within those sessions as they were by 2.5s GC perturbations. We suggest that such network-level effects on behavior reveal the widespread nature of taste processing, and the status of GC as just one (important) participatory node.
Despite affecting just one node of this large network of brain regions, our brief perturbations of GC uncover a temporally-specific role of cortex in driving (in a modulatory domain) orofacial behavior -a role that could not be discerned through wholesale disruption of activity. Our 0.5s perturbations reveal that GC contributes to the instigation of a gaping bout, but that GC plays no role in the maintenance of gaping once it begins. These data point at a dynamic flow of processing control within the larger taste network: modulatory signals propagate out of GC (signals that likely develop under the guidance of basolateral amygdala; Piette et al. (2012)) to influence the start of the motor program, which is then implemented and controlled by brainstem circuits. At its heart, the role of cortex in this model of taste processing has deep similarities to the role of neuromodulatory systems in the circuits underlying Aplysia feeding (Dacks and Weiss (2013)), leech swimming (Crisp and Mesce (2004)), control of gastric rhythms in the lobster and crab (Marder and Bucher (2007)), and rat whisking (Hattox et al. (2003)); in each, temporal aspects of rhythmic motor programs produced autonomously by a pattern generating circuit are influenced by descending signals.
The discreteness, coherence and inter-trial variability of GC ensemble activity states has attractor network-like features (Amit (1992), Hopfield (1982)). For one, attractor networks with multiple quasi-stable states can reproduce the sudden switches of activity seen in GC ensembles (Miller and Katz (2010)). In addition, the transition durations and state lifetime statistics of GC population dynamics are more in line with a dynamically switching attractor model than linear models of firing rate evolution (Jones et al. (2007), Sadacca et al. (2016)). Finally, nonlinear attractor-based circuits that exploit the noise inherent in neural processing more optimally perform the decision to ingest or expel a taste, which rats need no training to perform, than do linear integrating circuits (Miller and Katz (2013)). Our optogenetic protocol, with its mix of inhibitory and excitatory effects, presumably introduces a transient disruption in such attractor dynamics; such a perturbation is strong enough to transiently "knock" the network out of its quasi-stable sequence of states, but only if it hasn't already settled into the eventual, decision-related stable state.
The finding that the involvement of GC in the gape instigation process appears to last almost precisely 50ms past the calculated transition times can be attributed to one (or many) of several reasons. Firstly, transitions between quasi-stable states of GC processing, however discrete, are certainly not instantaneous -the time constants of neural firing ensure that there is some finite (albeit small) amount of time across which the ensemble makes the "jump" from one state of activity to another. In addtion, it is worth noting that both HMMs and change-point analysis techniques provide only a noisy estimate of state transition times, even if the transitions themselves were instantaneous. Finally, the change-point analysis model, by neglecting neural activity while the perturbation is happening, artifactually detects change-points close to the laser onset time on some of the trials, even if palatability firing actually began after the lasers were switched off. All of these factors likely play into the delay in gaping latency observed in the small subset of trials (28/168, 16%, Figure 7) with 0.7-1.2s perturbations where the lasers were switched on within 50ms of palatability change-point; while many of these trials reflect the non-instantaneous nature of transitions in ensemble activity patterns, (at least) some have artifactual change-points before laser onset without palatability firing actually emerging at that time.
In this study, we focused exclusively on gapes, the orofacial responses that rats make to expel aversive tastes from the oral cavity. Pilot attempts to implant EMG electrodes in deeper muscles that control licks (lateral tongue protrusions, LTPs) in response to palatable tastes resulted in unacceptable levels of distress for the animals. Although there are suggestions that gapes and LTPs can be produced by separate cortical mechanisms (Peng et al. (2015)), we consider that possibility unlikely for several reasons: 1) GC ensembles reflect the palatability of both appetitive and aversive tastes (Figure 3, Katz et al. (2001)) -even if palatability is modified by learning (Moran and Katz (2014)); 2) the latency and inter-trial variability of the onset of palatability-related ensemble activity is similar for palatable and aversive tastes (Sadacca et al. (2016)); and 3) there is considerable overlap in the brainstem circuits that underlie gapes and LTPs (Travers et al. (2000), Chen and Travers (2003), Venugopal et al. (2007), Moore et al. (2014)), resulting in similar latencies in the onset of LTPs and gapes after taste delivery (Travers and Norgren (1986)). These lines of evidence are consistent with the suggestion that cortex plays similar roles in the initiation of LTPs and gapes, which leads us to speculate that the transition of GC population activity to reflect stimulus palatability marks a shift in processing control, irrespective of the palatability of the tastant.
In summary, the balance of our results demonstrate a dynamic role for cortex in the processing of tastes; because this role involves ensemble activity states with variable trial-to-trial latencies, it cannot be discerned using standard analyses that average across trials. They reveal the importance of a unique moment in time that, despite being massively variable from trial to trial, denotes a reliable shift of processing control out of cortical circuits -a modulatory signal emerging (at least partly) from cortical circuits that is passed (presumably) to a brainstem central pattern generator. These results point at an attractor-like network of activity, potentially spread across interconnected brain regions, underlying the animal's decision to ingest or expel the tastant in the mouth -perturbations to this network can disrupt its functioning transiently, but only if it has not yet settled into the final, behaviorally-relevant stable state.  Figure 1: Experimental paradigm. A: Rats underwent surgeries for virus injections and implantation of opto-trodes and EMG electrodes. Post recovery from surgeries, they were given intra-oral infusions of Dil Suc (30mM Sucrose), Conc Suc (300mM Sucrose), Dil Qui (0.1mM Quinine-HCl) and Conc Qui (1mM Quinine-HCl) while GC neurons were inhibited by shining green (532nm) laser light on ArchT expressing GC neurons. B: Coronal slice from one of the animals in the study. ArchT expression (visualized by the GFP tag) is localized in GC. A lesion left by the tip of the opto-trode is visible in the middle of the GFP expressing region. C: Inhibition protocol used in the study. We delivered two types of optogenetic perturbations, short (0.5s) or long (2.5s), in separate experimental sessions.

A2
A3 A4 B1 B2 C D Figure 2: Impact of ArchT-mediated inhibition on GC neurons. A1-A4: Spiking rasters of an example unit in the 0.5s perturbation condition. Activity is robustly suppressed during laser stimulation. B1-B2: Spiking rasters of an example unit in the 2.5s perturbation condition showing clear inhibition during laser stimulation. C: Histogram of the change in firing rate (as a fraction of the rate on control trials) produced by 0.5s perturbations. The majority of neurons show robust suppression when perturbed; a small group of neurons increase their firing rate in response to perturbation, presumably due to network-level effects. The largest fraction of suppressed neurons are seen when perturbation is delivered in the middle of the GC taste response, from 0.7s to 1.2s post stimulus. D: Histogram of the change in firing rate (as a fraction of the rate on control trials) produced by 2.5s perturbation. Almost all neurons are affected by the perturbation: the large majority are suppressed, but a small minority show elevated firing rates in response to perturbation. Figure 3: Impact of optogenetic perturbations on palatability relatedness of the firing of GC neurons. A: Coefficients obtained from the regression of trial-averaged firing rates on palatability ranks of the taste stimuli. The solid lines depict the mean regression coefficient across time while coefficients significantly different from 0 at the 5% level are marked by dots. The dashed lines are logistic sigmoid fits to the mean regression coefficient in each condition. Disruption of GC firing for 2.5s wipes out the entirety of the palatability response, which never differs significantly from 0. B: The post stimulus latency (blue bars) and magnitude (red bars) of the peak (95% of the asymptote) of the sigmoid fits in A. Error bars denote 95% Bayesian credible intervals and indicate statistical significance at the 5% level if they are not overlapping. On control (laser off) trials, GC neurons asymptote to peak palatability firing ∼ 0.8s post stimulus. The 2.5s perturbation, by disrupting the palatability response completely, is fit by a flat sigmoid whose peak magnitude overlaps 0. C: The post stimulus latency (blue bars) and magnitude (red bars) of the peak (95% of the asymptote) of the sigmoid fits in the 0.5s perturbation condition. Error bars denote 95% Bayesian credible intervals and indicate statistical significance at the 5% level if they are not overlapping. On laser off trials, GC representation of palatability peaks ∼ 0.8s after taste delivery, identical to the 2.5s perturbation control trials in B. Perturbations early (0-0.5s) and in the middle of the taste response (0.7-1.2s) delay the peak of palatability firing by ∼ 0.5s; the magnitude of this peak, however, is the smallest for the middle perturbation. Perturbations late in the taste trial (1.4-1.9s), after palatability firing is mostly over, expectedly have no impact and align with the control trials. Probability of gaping is calculated as the total posterior probability of 4-6Hz movements. A1-A2: Two representative Conc Qui trials. The animal's mouth movements can be seen as bursts of EMG activity (blue) following taste delivery -the onset of gaping, as detected on video, is marked. The envelope of the EMG signal (black line) was subjected to spectrum analysis. B1-B2: Two representative Conc Qui trials. The probability of gaping from BSA (black line) matches up with individual gapes (grey notches) picked by a previously published quadratic classifier that achieved 75% accuracy. C1-C2: BSA (solid line) and the quadratic classifier (dotted line) produce similar estimates of trial-averaged probability of gaping in response to Dil Qui (C1) and Conc Qui (C2) on a set of control (laser off) trials. D: The probability of gaping from BSA rises reliably just before the first gape. Gaping probability was averaged across trials aligned by the time of the first gape, detected either on video (black) or by the quadratic classifier (grey). The black dashed line (0 on the x axis) indicates the occurrence of the first gape. E: KL divergence between the probability of gaping to Conc and Dil Qui (higher values indicate larger differences in their gaping distributions, same trials as in B). As expected, the distributions of gaping probability on Conc and Dil Qui trials are initially similar (while nonspecific investigative licks happen) and diverge out at ∼1s post stimulus once gaping begins. F: The cumulative sum of the KL divergence in E across time. The jump in KL divergence around the mean onset time of gaping is seen as a change in slope of its cumulative sum. We fit two straight lines to the cumulative sum and pick their intersection as the mean onset of gaping across this set of trials. Figure 5: Onset times of 4-6Hz aversive orofacial behaviors (gapes). The bars depict the mean gape onset times while the extent of their 95% Bayesian credible intervals are shown by the error bars. Non-overlapping error bars depict statistical significance at the 5% level. A: Onset of aversive orofacial behaviors in the control trials. The 2.5s controls show a delayed onset, likely due to spillover effects from the long optogenetic perturbation. B: Onset of aversive orofacial behaviors in the experimental trials. Perturbations early (0-0.5s) and in the middle (0.7-1.2s) of the taste response bring about a similar delay in the onset of gaping as the 2.5s perturbation. Gaping is unaffected if GC neurons are disrupted late in the trial (1.4-1.9s) The transitions cannot happen during the perturbation, from 0.7s to 1.2s. A2: Distribution of time of transition to the palatability-related ensemble state on Conc Qui trials. The transitions cannot happen during the perturbation, from 0.7s to 1.2s. B: Coefficients from the regression of trial-averaged firing rates of GC neurons on palatability ranks of the taste stimuli. The lines depict the mean regression coefficient across time while coefficients significantly different from 0 at the 5% level are marked by dots. The coefficients are significantly different from 0 only on the trials where the palatability-related ensemble state appeared before the perturbation started. C: Onset of aversive orofacial behavior when GC is perturbed from 0.7s to 1.2s post stimulus. The onset of gaping is delayed if the perturbation begins before palatability information has appeared in ensemble activity. A subset of trials, with palatability transition times in the 50ms prior to 0.7s, appear to produce a delay in gaping behavior despite the emergence of the ensemble state before the start of the optogenetic perturbation. Dropping these trials from the analysis revealed the gaping happens earlier than on control trials if GC ensembles shift to palatability firing before the lasers are switched on.