## New Research In

### Physical Sciences

### Social Sciences

#### Featured Portals

#### Articles by Topic

### Biological Sciences

#### Featured Portals

#### Articles by Topic

- Agricultural Sciences
- Anthropology
- Applied Biological Sciences
- Biochemistry
- Biophysics and Computational Biology
- Cell Biology
- Developmental Biology
- Ecology
- Environmental Sciences
- Evolution
- Genetics
- Immunology and Inflammation
- Medical Sciences
- Microbiology
- Neuroscience
- Pharmacology
- Physiology
- Plant Biology
- Population Biology
- Psychological and Cognitive Sciences
- Sustainability Science
- Systems Biology

# Modeling the effect of locus coeruleus firing on cortical state dynamics and single-trial sensory processing

Contributed by Nikos K. Logothetis, August 20, 2015 (sent for review June 22, 2015)

## Significance

Understanding what makes a cortical neuron fire is fundamental to understand how neural circuits operate. We used simultaneous recordings from the locus coeruleus (the neuromodulatory nucleus releasing norepinephrine) and from the somatosensory cortex to formulate a mathematical model explaining how cortical responses originate from the interplay of the sensory drive that cortical neurons receive, the spontaneous dynamics of cortex, and the effect of neuromodulation. Our work provides a hypotheses about how the temporal structure of locus coeruleus burst firing regulates the amplitude and timing of changes in cortical excitability and may selectively amplify responses to salient sensory stimuli. It also suggests that downstream circuits may better decode the activity of a cortical sensory network after estimating its neuromodulatory state.

## Abstract

Neuronal responses to sensory stimuli are not only driven by feedforward sensory pathways but also depend upon intrinsic factors (collectively known as the network state) that include ongoing spontaneous activity and neuromodulation. To understand how these factors together regulate cortical dynamics, we recorded simultaneously spontaneous and somatosensory-evoked multiunit activity from primary somatosensory cortex and from the locus coeruleus (LC) (the neuromodulatory nucleus releasing norepinephrine) in urethane-anesthetized rats. We found that bursts of ipsilateral-LC firing preceded by few tens of milliseconds increases of cortical excitability, and that the 1- to 10-Hz rhythmicity of LC discharge appeared to increase the power of delta-band (1–4 Hz) cortical synchronization. To investigate quantitatively how LC firing might causally influence spontaneous and stimulus-driven cortical dynamics, we then constructed and fitted to these data a model describing the dynamical interaction of stimulus drive, ongoing synchronized cortical activity, and noradrenergic neuromodulation. The model proposes a coupling between LC and cortex that can amplify delta-range cortical fluctuations, and shows how suitably timed phasic LC bursts can lead to enhanced cortical responses to weaker stimuli and increased temporal precision of cortical stimulus-evoked responses. Thus, the temporal structure of noradrenergic modulation may selectively and dynamically enhance or attenuate cortical responses to stimuli. Finally, using the model prediction of single-trial cortical stimulus-evoked responses to discount single-trial state-dependent variability increased by ∼70% the sensory information extracted from cortical responses. This suggests that downstream circuits may extract information more effectively after estimating the state of the circuit transmitting the sensory message.

Responsiveness of cortical sensory neurons is state dependent. In other words, the neural responses to a sensory stimulus do not only depend on the features of extrinsic sensory inputs but also on intrinsic network variables that can be collectively defined as the “network state” (1). Cortical sensory neurons receive information about the external world from peripheral receptors via feedforward sensory pathways. However, the abundance of recurrent and feedback connectivity (2) may generate ongoing activity that shapes the background on which the afferent information is processed (3). In addition, neuromodulatory inputs from neurochemically specialized brain nuclei that are not part of the direct spino-thalamo-cortical pathway can modulate the dynamics of cortical networks (4), as well as control the animal’s behavioral state. The concurrent integration of information about the external world and about internal states is likely to be central for computational operations of cortical circuits and for the production of complex behavior, yet its mechanisms and implications for neural information processing are still poorly understood.

The locus coeruleus (LC) is a brainstem neuromodulatory nucleus that likely plays a prominent role in shaping cortical states via a highly distributed noradrenaline release in the forebrain (5). In particular, the LC contributes to regulation of arousal and sleep; it is involved in cognitive functions such as vigilance, attention, and selective sensory processing (5⇓–7); and it modulates cortical sensory responses and cortical excitability (8).

Here, we investigated how LC firing influences spontaneous and stimulus evoked cortical activity by performing simultaneous extracellular recordings of spontaneous and somatosensory-stimulation–evoked neural activity in the primary somatosensory cortex (S1) and in both ipsilateral LC (i-LC) and contralateral LC (c-LC) in urethane-anesthetized rats. We first use these data to investigate the statistical relationships between the temporal structure of LC firing and the changes in cortical excitability, and we then construct a dynamical systems model of the temporal variations of cortical multiunit activity (MUA) that describes quantitatively the dynamic relationships between LC firing, spontaneous cortical activity dynamics, and cortical sensory-evoked responses. We use this model to study how a specific neuromodulatory input may influence the information content and the readout of cortical information representations of sensory stimuli.

## Results

### Neurophysiological Recordings from Somatosensory Cortex and LC in Urethane-Anesthetized Rats.

To investigate the relationship between the dynamics of noradrenergic cell firing and of ongoing and stimulus-evoked cortical activity, we recorded simultaneously the multiunit spiking activity in S1 and in LC of urethane-anesthetized rats (*n* = 4). MUA in S1 was recorded using a linear electrode array (28 contacts with 50-μm spacing) covering most of the cortical depth and placed in the S1 hindpaw representation in the hemisphere contralateral to the stimulated paw. Hereafter, we analyzed the mass cortical MUA obtained pooling spiking activity from all cortical electrodes (abbreviated as “cortical MUA”). The LC MUA was recorded using two single electrodes each placed in the LC core, i-LC and c-LC to the S1 recording site. We recorded data both during stretches of spontaneous neural activity and during intermittent somatosensory stimulation. Specifically, trains of five electrical pulses were delivered to the hindpaw; the pulses were presented at three different frequencies (15, 25, and 35 Hz) and three different amplitudes (1, 3, and 5 mA), and each of the nine stimulus types were delivered in random order with a 10-s interval and were repeated 20–30 times per session.

In absence of sensory stimulation (spontaneous periods), cortical spiking activity was highly synchronized across neurons and underwent slow variations in excitability, with typically one to four peaks of spontaneous firing per second alternating with periods of reduced firing (Fig. 1 *A* and *B* and Fig. S1 *A* and *B*). Consequently, the power spectrum of spontaneous cortical MUA (Fig. 1*C*) had the highest power at low frequencies. Spontaneous MUA of i-LC and c-LC was characterized by phasic elevations of firing (called “bursts” hereafter) alternating with periods of reduced activity (Fig. 1 *A* and *B*). The bursts of LC MUA happened typically four to seven times per second, and thus were more frequent than the peaks of spontaneous cortical firing. Power spectra of LC MUA (Fig. 1*C*) showed a local peak in the 4- to 7-Hz range, capturing the rhythmic nature of the timing of LC bursts. S1 neurons responded to the pulses of somatosensory stimulation (Fig. 1*B* and Fig. S1 *C* and *D*) with a short-latency discharge (likely reflecting primarily spinothalamic input) and a longer-latency sustained component (likely reflecting also nonspecific input and engagement of recurrent or top-down networks). After somatosensory stimulation, both LCs fired before the S1 longer-latency response (Fig. 1*B* and Fig. S1 *C* and *D*).

### How the Dynamics of LC Firing Influences Cortical State Changes.

To understand the role of noradrenergic cell firing in shaping cortical dynamics, we first investigated whether there were systematic temporal relationships between them during spontaneous activity. (Unless otherwise stated we hereafter report mean ± SEM as statistics computed over the dataset.)

We computed the cross-correlogram between the time course of i-LC firing [which was the LC that—in agreement with anatomical connectivity (9) and other results reported in *Supporting Information*—influenced ipsilateral cortical MUA the most] and that of cortical MUA, finding that it peaked when i-LC activity was shifted back by ∼20–25 ms (Fig. S2*A*). In addition, we computed—using transfer entropy (TE) (see ref. 10 and *Supporting Information*)—the amount of causation exerted by i-LC MUA on cortical MUA, and we compared it with the amount of causation exerted by cortex on i-LC. We found that the TE from i-LC to cortex was significant (*P* < 0.001, *t* test) and was 2.9 ± 1.1 times larger (*P* < 0.001; paired *t* test) than in the reverse direction (Fig. 1*D*; see also Fig. S2). Together, these results support the hypothesis that LC firing exerts a causal influence on cortex within few tens of milliseconds.

To further investigate how the temporal structure of LC firing may influence that of cortical activity, we then computed the correlation between the frequency spectrum of i-LC MUA and that of S1 MUA (both computed from the same 1.5-s spontaneous activity window). To focus on genuine time variations of LC firing, we normalized the power of i-LC MUA to its time average in the trial. We found (Fig. 1*F*) a significant (*P* < 0.01, paired *t* test) positive correlation between the low-frequency (1–10 Hz) power of i-LC MUA and the power of delta-band (1–4 Hz) cortical MUA. This positive correlation suggests that one effect of the low-frequency rhythmicity of LC is to increase the power of low-frequency synchronized fluctuations of cortical excitability. In other words, the low-frequency rhythmicity of LC MUA amplifies low-frequency cortical state variations. We note that the delta cortical power did not correlate (*P* > 0.2, *t* test) with LC averaged firing rate (Fig. S3*A*), and thus the amplification of cortical slow oscillations was specifically related to low-frequency rhythmicity of LC activity. However, there was a significant (*P* < 0.01, paired *t* test) negative correlation between LC firing rate and S1 MUA power in the 5- to 20-Hz range, suggesting that the level of tonic LC firing varies across cortical states (11).

Finally, we investigated the relationship between the timing of i-LC firing and of the cortical delta-band state fluctuations. We expressed the timing of low-frequency fluctuations between periods of low and high cortical excitability using the phase of delta-band cortical MUA, with phase π and 0 representing the peak and trough of cortical firing, respectively (Fig. 1*E*). We found that—consistent with results obtained during NREM sleep (8)—i-LC was significantly phase-locked (Rayleigh test, *P* < 0.05) to the cortical delta rhythm, and that it fired preferentially one-eighth of cycle (π/4 radians) ahead of peak cortical excitability. This suggests that a burst of LC firing (such as, e.g., those plotted in Fig. 1 and Fig. S1) facilitates the transition to a state of higher cortical excitability after a few tens of milliseconds.

### A Dynamical Systems Model of the Relationship Between LC Firing and Cortical State Changes.

Our goal was to characterize quantitatively if and how the dynamics of the spontaneous low-frequency changes in cortical excitability interacts with the firing of LC neurons and with the afferent sensory drive to generate state-dependent cortical representations of stimuli.

Extending the work of ref. 12, we described quantitatively the time course of cortical MUA with a simple discrete-time dynamical system inspired by the FitzHugh–Nagumo model (13). In brief, the model (schematized in Fig. 2*A* and Fig. S4, and denoted hereafter as the “true-LC model”) describes cortical activity in an effective (rather than biophysically detailed) way. It contains activity variables that can be interpreted as summarizing the total rate of excitatory and inhibitory neurons in a cortical network and it contains terms of self-excitation and inhibition that lead to ongoing intrinsic fluctuations of the activity, as well as an input term describing the stimulus drive to the cortical network following the application of a somatosensory stimulus. Importantly, it also models the potential causal effect of LC on cortex by including as input the MUA time series from both i- and c-LC (shifted in the past by a lag of 20 ms, whose value optimized model performance; *Supporting Information*). The best-fit parameters define our estimate of the cortical state and are used to predict the cortical activity afterward. This model contains both “activity state” model variables (such as S1 and LC firing) that describe cortical dynamics persisting on shorter timescales of few tens to few hundred milliseconds, and “dynamic-state” variables (such as the couplings in the dynamical systems model that specify how cortical MUA evolve in time) that describe properties of network dynamics on a timescale of seconds (12).

To determine quantitatively cortical states, we first fitted the model to periods of spontaneous activity by minimizing (over a 1.5-s “state detection” or “fit” period) the error between the S1 MUA and its model prediction. To test how well the model kept tracking cortical dynamics, we then generated a prediction of cortical MUA in a postfit period (by solving the dynamics of the model with the best-fit parameters and the initial conditions given by the activity at the end of the fit period) for a further “prediction period” (until 500 ms after the end-of-fit period when predicting continuation of spontaneous activity, and until 300 ms after the end of paw stimulation when predicting responses to stimulation). To evaluate the specific effect of the time course of LC firing on the dynamics of synchronized cortical states, we compared the performance of the true-LC model (in which the time evolution of cortical MUA depended both on past cortical activity and past LC activity) to that of a second model, named “no-LC” model and similar to the one introduced in ref. 12, in which evolution of cortical MUA depended only on its own past but not on LC firing (*Supporting Information*). Examples of fits to individual periods of cortical spontaneous and stimulus-evoked activity are reported in Fig. 1 *A* and *B* and Fig. S1 *A–D*. Although both models followed reasonably well cortical firing dynamics, the true-LC model fitted and predicted it better. In particular, the true-LC model generated low-frequency fluctuations of cortical MUA with higher power than the no-LC model, and this led to a better match with the real cortical MUA (Fig. 1 *A* and *B* and Fig. S1 *A–D*). The true-LC model predicted better also the single-trial cortical response to the somatosensory stimulus (Fig. 1*B* and Fig. S1 *C* and *D*), suggesting that its more accurate representation of cortical states leads to a better prediction of the interaction between stimulus drive and ongoing activity that may generate the single-trial response to the stimulus.

We quantified systematically the ability of the different models to represent the dynamics of ongoing spontaneous cortical activity. We found (Fig. 2*B*) a significant (*P* < 0.001, paired *t* test) advantage (35% reduction of fit error) of the true-LC model over the no-LC model in describing cortical dynamics during the fit period. This advantage was genuine and not simply due to the larger number of parameters in the LC model (*Supporting Information*). Moreover, the true-LC model predicted better the MUA in postfit periods of spontaneous activity, with an increase (*P* < 0.01; paired *t* test) of ∼90% of the Pearson correlation between real and model-predicted cortical MUA (Fig. S5*A*). The true-LC model also predicted better than the no-LC model the trial-averaged amplitude of the cortical MUA evoked by each stimulus (signal correlation between model and data; Fig. S5*B*; *P* < 0.001, paired *t* test) and the trial-to-trial variability of single-trial responses around their mean (noise correlations between model and data; Fig. 2*C*; *P* < 0.01, paired *t* test). In sum, the true-LC model predicted better both spontaneous and stimulus-evoked S1 dynamics. We performed further analyses of our models (reported in full in *Supporting Information*) that established that the prediction of the variations across trials of cortical stimulus-evoked MUA obtained from both the no-LC and the true-LC model was much better than the prediction obtained from simpler proxies of cortical states such as the synchronization index (S.I.; Fig. S6*C*) or the cortical delta-phase or MUA amplitude at the time of stimulus application (Fig. 2*C*), demonstrating the value of dynamical systems in describing cortical states.

To characterize the aspects of cortical dynamics that were better described by the true-LC model, we calculated the difference between the power of real cortical activity and that of the model at each frequency during spontaneous cortical dynamics (normalized dividing it by the total power in the 1- to 10-Hz range). This difference was positive for both models but was much larger for the no-LC model, especially in the delta frequency band (Fig. 3*A*), suggesting that the no-LC model generates too little delta-band power in comparison with real cortical MUA, and that introducing LC firing into the model increases the power of the low-frequency cortical fluctuations that it can generate, thereby improving the accuracy of the description of cortical dynamics. We then computed—as we did in real data—the correlation between spectra of i-LC MUA and of model S1 MUA. We found (Fig. S7*C*) a positive (*P* < 0.01, *t* test) correlation between the power of the i-LC MUA at low frequencies (1–10 Hz) and the power of model delta-band cortical MUA, as in real data. This suggests that the model captures the key role of low-frequency LC rhythmicity in increasing the power of low-frequency synchronized cortical state variations (see *Supporting Information* for further analyses). The reason was that the best-fit parameters (Fig. S8) of the nonlinear coupling term between i-LC and cortical MUA had a sign that allowed LC bursts to either increase or decrease cortical excitability depending on whether cortical MUA is in the ascending or descending phase of its up-down state dynamics, thus adding oscillatory power to cortex (Fig. S9).

Finally, we found that the phase coherence between real and model cortical MUA was higher (*P* < 0.001, paired *t* test) for the true-LC than for the no-LC model (Fig. 3*B*), because the true-LC model uses the timing of LC bursts to better predict the timing of peaks and troughs of cortical excitability.

### Extracting Sensory Information from State-Dependent Cortical Codes: Discounting State-Dependent Trial-to-Trial Variability.

The large trial-to-trial variability of cortical responses is a major obstacle in the readout of neural codes (14). Our results above—and several previous studies—show that variations in cortical state (12, 15⇓–17) and neuromodulation (4) contribute to trial-to-trial cortical response variability. Here, we explore the consequences of successfully modeling the neuromodulatory cortical state for the decoding of cortical responses to sensory stimuli.

Given that our models can partly account for and predict the trial-to-trial variations of stimulus-evoked MUA responses, these model predictions can be used to increase the information about stimuli extracted from cortical activity. More sensory information can be obtained by simply subtracting out the estimated state-induced trial-to-trial variations of these responses. This idea is illustrated with real data in Fig. 4*A*, which shows the time course of the trial-averaged stimulus-evoked S1 MUA to two stimuli of different intensity in a single example trial. This trial elicited a firing rate whose strength was in-between the mean rates of these two stimuli. This intermediate- strength response could conceivably have arisen either in response to the weakest stimulus when the network was in excitable state or in response from the strongest stimulus when the network was less excitable. This ambiguity can be resolved by computing the trial-to-trial variability predicted by the network state. In this example, the predicted variability was positive (indicating that the network was in a more excitable state). The subtraction of the predicted variability from the single-trial response (black downward arrows in Fig. 4*A*) produces a “variability-discounted” response much closer to the trial-averaged response of stimulus presented in that trial (and thus much easier to decode) than the original response.

To quantify the effect of the reduction of variability on stimulus discrimination, we computed (see ref. 14 and *Supporting Information*) the information about which somatosensory stimulus was presented that is available in the original cortical MUA, and the information that can be extracted when subtracting out the state-dependent variability. Note that the prediction of state-dependent variability was obtained by the model without any knowledge about which stimulus was presented, and that subtracting variability may not always be the optimal way to discount it. Thus, the difference between these two information values is a robust lower bound to the information that can be gained by discounting the state-dependent variability captured by the models.

To avoid confounds resulting from differences in the pulse application times, we computed only information between pairs of stimuli with equal pulse frequency, but different intensity of the delivered electrical current. The time course of the information instantaneously available in these neural activities after the end of stimulus, averaged over all pairs of stimuli and sessions, is reported in Fig. 4*B*. We found that discounting the variability of MUA with the output of both true-LC and no-LC models leads to obtain significantly (*P* < 0.001, paired *t* test) more information than with the nondiscounted original MUA. The gain in information due to state variability discount was particularly large in the late (125–200 ms) sensory response component, which has been shown to be particularly relevant for behavior (18). Across the 125- to 200-ms late-response period, the gain of information of the denoised responses over the information present in the original MUA was in the range of 32–71% for the true-LC model and 18–39% for the no-LC model. The advantage of the true-LC model over the no-LC model was statistically significant (*P* < 0.05, paired *t* test). Importantly, the information gained by discounting the variability estimated by the dynamical state models was much larger than the information that would be obtained by simply subtracting out the MUA at the time of stimulus application of by predicting state with S.I. or the cortical delta-phase or MUA amplitude at the time of stimulus presentation (Fig. S10*A*). This shows the worth of dynamical systems for extracting information from neural data.

### Modeling the Effect of LC Phasic Activation on Cortical Responses to Stimuli.

We finally used our model to study how phasic activation of LC may shape the cortical information representation of sensory stimuli. Phasic activation of LC refers to brief transient bursts of firing of LC neurons, which may happen spontaneously or be brought about by electrical stimulation of LC or by a salient external event (5).

We created a series of simulated single-trial cortical responses to somatosensory stimuli using the true-LC model with the best-fit parameters and stimulus-drive functions taken from all trials in real data. For all such simulated trials, we ran two sets of calculations. In the first set, describing normal spontaneously occurring states, we ran the model using the true-LC time course in that trial as neuromodulatory input. In the second set, we simulated the effect of phasic activation of LC by adding to the true-LC dynamics a brief burst of firing (Fig. 5*A*) whose amplitude and timing was chosen to maximize the amplitude of the resulting stimulus-evoked response. We then inferred the specific effect of phasic LC activation by comparing the stimulus-evoked cortical MUA across these two simulations. We found (Fig. 5*B*) that phasic LC activation increased the time-averaged amplitude of the MUA in the 125- to 200-ms post–end-of-stimulus late component that gives maximal stimulus information. We also observed a reduction of trial-to trial variability in response peak time and in rise time during the late response window when using phasic LC activation before stimulus presentation (Fig. 5*C*). Significant (*P* < 0.001, paired *t* test) reduction of variability was obtained even in the presence of trial-to-trial jitter between time of phasic LC activation and of stimulus application, as long as the jitter range was below 110 ms for peak time and below 40 ms for rise time (Fig. S11).

The increased response to weaker stimuli and the increased precision of cortical spike times following phasic electrical stimulation of LC has been empirically reported (6, 19). Our model provides a possible mechanistic explanation of such previous findings. Our model predicts that, in a nonvigilant state (as is the rat state in our study), cortex will be more consistently excitable if the stimulus is presented shortly after the additional phasic LC activation, thereby making cortex more responsive particularly to weaker stimuli and quenching part of the fluctuations of excitability that cause variability of stimulus-evoked responses.

## Discussion

Recently, much research has focused on how brain states affect the dynamics of cortical sensory responses (1, 15, 16, 20). In particular, statistical models have been developed to describe the relationships between indicators of network states and the spike counts observed in response to a stimulus (15⇓–17). Our study adds critical insights about the dynamic interplay of spontaneous and stimulus-evoked activity. By formalizing state dependence using a dynamical system (12), we were able to capture important temporal relationships—for example, the relationship between periodic burst of LC firing and the amplitude and timing of changes in cortical excitability—which are not captured by spike count models that average the neural response over time. Moreover, by outperforming a widely used frequency analysis-based predictor of behavioral responses and of cortical responsiveness from neural recordings—the phase of low-frequency oscillations (21⇓–23)—the dynamical systems approach shows the ability to capture complex nonperiodic temporal features that are relevant for sensory processing and that cannot be fully described by a frequency-band analysis. This underscores the promise of the dynamical systems approach to elucidate how neural circuits shape sensory processing.

An important addition of our work to previous models of state dependence was the inclusion of the contribution of an important neuromodulator—the noradrenergic system. Our results support the hypothesis that the temporal structure of LC firing causally influences cortical dynamics and that, in particular, LC bursts may play a key role in regulating both the amplitude and timing of low-frequency cortical oscillations. Our model provides two predictions that can be tested by future interventional experiments designed to specifically manipulate the temporal structure of LC bursts: that their low-frequency rhythmicity increases the amplitude of cortical delta fluctuations and that their timing regulates when cortex is maximally responsive to stimuli and influences the spike timing reliability of cortical neurons.

Although slow oscillations have a neocortical origin (24), our results raise the possibility that LC may be a key component of a larger network including thalamocortical loops (25) that may regulate slow oscillation and mediate their large-scale synchronization. LC phasic activity may modulate slow oscillations via direct coerulear-cortical projections or indirectly via temporal synchronization of activity within multiple members of the ascending system or cortico-thalamic loops. The role of LC in this potentially wide network could be investigated by generalizing this work to model the interaction between LC and multiple simultaneously recorded cortical and thalamic regions, and the resulting model predictions could be directly tested by future interventional experiments.

Our data were obtained in anesthetized animals. However, this still allows speculating about the influence of phasic LC firing when salient stimuli are presented during a nonvigilant state. Phasic LC activity may selectively facilitate responses to task-relevant processes and filter away responses to irrelevant events (7) or orient attention to salient events (5). Fluctuations in cortical excitability may be a key mechanism for the selection of salient variables (23) as they permit selective amplification or attenuation of stimulus responses depending on whether stimuli are presented at more or less excitable states. Thus, the ability of the low-frequency rhythmicity of LC bursts to enhance the power of fluctuations in cortical excitability may enhance salient variables. Our work highlights the importance of timing of LC burst: suitably timed LC burst (for example, triggered by an alerting stimulus) can very rapidly trigger transitions into excitable cortical states, which in turn decrease the threshold for cortical responses and thus dynamically facilitate the processing of salient or attended events. That LC bursts may modulate dynamically the threshold for cortical processing and regulate stimulus selection was proposed in connectionist models (26). Our work identifies potential neural mechanisms and reveals the precise timescales of this phenomenon.

The state dependence of neural responses profoundly constrains how neural circuits exchange information. State dependence may either force neurons to transmit information only using codes that are robust to state fluctuations (e.g., relative firing rates), or may force downstream neurons to gain information about the state of the networks sending the sensory messages and then to use the knowledge of state to properly interpret neural responses. Our results suggest that the latter information transmission scheme is feasible, because detecting state by either monitoring the dynamics of cortical ongoing activity alone or by also monitoring the dynamics of noradrenergic modulation substantially increased the amount of information about sensory stimuli in the late response components relevant for behavior. Given that ongoing dynamics appears correlated over relatively large spatial scales and that norepinephrine is released diffusely, it seems conceivable that state variables are shared by other neurons in the same or a nearby region. Because most cortico-cortical synapses connect proximal neurons (27), the state variables will often be shared by the neuron sending a message and the neuron receiving it; thus, state-dependent decoding may work well for many cortical computations.

## Materials and Methods

These data were recorded from four male Sprague–Dawley rats (250–350 g). All experiments were approved by local authorities (Regierungspräsidium Tübingen, Germany, Referat 35, Veterinärwesen), in accordance with the regional animal welfare committee pursuant to §15 of the German Animal Welfare Act (Kommission nach §15 des Tierschutzgesetzes), and fully complied with Directive 2010/63/EU of the European Union and of the council on the protection of animals used for scientific purposes. MUA was simultaneously collected from 28 electrodes (array spacing, 50 μm; impedance, 200 kΩ) placed in the paw representation area of the right-hemisphere S1 cortex (S1) while the rats were anesthetized with urethane. At the same time, MUA was also recorded from both left and right LC using single microelectrodes. We recorded responses to nine different trains of electrical left hindpaw stimulations (three different frequencies and three different amplitudes) and during spontaneous activity before the stimulus application.

A finite-difference dynamical systems model of cortical S1 MUA (Eqs. **S1**–**S5**) was fit to the data by choosing model parameters minimizing the fit error (Eq. **S9**) between real S1 MUA and its model prediction during a 1.5-s fit period of spontaneous activity. The model was also used to predict the time course of cortical MUA for further periods of spontaneous or stimulus-evoked activity. We described the data using two models: the true-LC model that contains a coupling term between LC activity and cortical dynamics (Eq. **S5**), and a no-LC model that does not have this coupling term. Information and TE were computed from Eqs. **S11**–**S13** using a brute-force estimation of the neural response probability discretized into a small number of bins and using bias corrections (28).

Full details on neurophysiological procedures and analysis are described in *SI Materials and Methods*.

## SI Materials and Methods

### Ethical Statement.

The data used in these analyses were recorded in full accordance with both local and European Union regulations. The data from four male Sprague–Dawley rats (250–350 g) were used in this study. The animals were housed on a 12-h light/dark illumination cycle with constant access to food and water. All experimental procedures were approved by the local authorities (Regierungspräsidium) and were in full compliance with the European Parliament and Council Directive 2010/63/EU on the protection of animals used for experimental and other scientific purposes.

### Brief Overview of Neurophysiological Recordings from Somatosensory Cortex and Locus Coeruleus in Anesthetized Rats.

Multiunit spiking activity was simultaneously collected from 28 electrodes (array spacing, 50 μm; impedance, 200 kΩ) placed in the paw representation area of the right-hemisphere primary somatosensory cortex (S1) while the animals were anesthetized with urethane. At the same time, multiunit spiking activity was also recorded from both the left and right locus coeruleus (LC). We recorded responses to nine different trains of electrical paw stimulations (three different frequencies and three different amplitudes, 15–30 trials per stimulus). Importantly, spontaneous activity in the few seconds preceding the stimulus presentation was also collected.

### Surgery and Electrophysiological Recording.

Rats were deeply anesthetized with urethane (1.5 g/kg, i.p.) and the surgery begun after the lack of withdrawal reflex in response to manual noxious pinch on the paw. The heart rate and blood oxygenation were monitored using pulse oxymeter (NoninR 8600V; Nonin Medical, Inc.), and oxygen was provided to maintain the blood oxygenation level above 80%. The body temperature was maintained at 37 °C throughout the experiment by a rectal probe and a feedback-controlled thermal blanket. The skull was exposed and anterior-posterior angle between the bregma and the lambda skull marks was set at 0°. The final recording coordinate in the primary somatosensory cortex (S1HL) was determined by mapping the highest epidural evoked potential to a single biphasic (0.5 ms, 5 mA) pulse delivered to the contralateral hindpaw by a pair of 26-gauge needles inserted s.c. [the typical coordinate was as follows: anteroposterior (AP) = −1.0 mm from bregma; mediolateral (ML) = 3 mm]. An array with 32 channels and 50-μm interchannel spacing (Neuronexus Technologies) was inserted in the mapped site approximate 1.8 mm below the dura matter, ensuring that the upper and bottom channels in the array displayed the inverse polarity. Bilateral recordings of LC were performed using extremely thin (5 μm) microelectrodes (Carbostar-3; Kation Scientific). The electrodes targeting the LC ipsilateral and contralateral to the recorded S1HL entered the brain with mirrored angles of 15° (AP = −4.0 mm from lambda; ML = 1.2 mm and 6.0 mm from the dura surface) and −15° (AP = −0.5 mm from lambda; ML = 1.2 mm and 6.5 mm from the dura surface) to the normal of the anterior-posterior axis, respectively. The final electrode coordinate in each LC was electrophysiologically guided by distinct activity pattern of norepinephrine (NE) neurons (29) and the NE nature of the recorded multiunit activity (MUA) was confirmed by i.p. injection of clonidine at the end of the experiments.

The extracellular signals were amplified (5,000×) using a custom-made 32-channel preamplifier and the Alpha Omega multichannel processor (MPC Plus; Alpha Omega Co.). The signals were digitized at 24 kHz using CED Power1401mkII converter and Spike2 data acquisition software (Cambridge Electronic Design).

### Somatosensory Stimulation.

For somatosensory stimulation, we applied electrical footshocks via two stainless-steel needles placed s.c. ∼1 cm apart in the paw of the hindlimb contralateral to the recording sites. Electrical current was delivered using a biphasic stimulus isolator (BSI-1; Bak Electronics). The stimulation parameters were digitally controlled by the Spike2 software (Cambridge Electronic Design) and transmitted to the current source via digital-to-analog converter built-in to the data acquisition unit CED Power 1401mkII (Cambridge Electronic Design). The voltage, which is necessary to pass the current between the needle tips, was monitored on the oscilloscope via a custom-built voltage output unit and compared with the digital input. The stimulation consisted of trains of five rectangular biphasic pulses (0.5-ms pulse duration) delivered to the hindpaw; the pulses were presented at three different frequencies (15, 25, and 35 Hz) and with amplitude of either 1, 3, or 5 mA. Each of these possible nine stimulus types (each made by a combination of a given frequency and a given amplitude) were applied in random order and repeated 20–30 times per session with a 10 ± 1-s jittered interval between stimuli.

### Perfusion and Histology.

At the end of experiment, the rat was euthanized with a lethal dose of sodium pentobarbital (100 mg/kg, i.p.; Narcoren; Merial GmbH) and perfused transcardially. Serial 60-μm-thick coronal sections were cut on a horizontal freezing microtome (Microm HM 440E), Nissl-stained, and examined under microscope.

### MUA Extraction.

For S1 channels, spike times were first detected from the bandpassed (400–3,000 Hz, fourth-order Butterworth Filter) extracellular potential in each electrode by threshold crossing (>3 SD). Mass MUA of S1 cortex was obtained by pooling together the trains of all recorded spikes form all electrodes. From the mass cortical MUA time series, we obtained a smoothed mass MUA times series, denoted as *i*th time bin was a weighted average of the MUA in the previous 100 time bins). The time series of the smoothed mass cortical MUA was the quantity that we fitted to the dynamical systems models of cortical state changes (see below).

Given that we placed only one extracellular electrode in each LC, the computation of MUA

### The Dynamical Systems FitzHugh–Nagumo Model of the Interaction Between Ongoing Cortical State Dynamics, Noradrenergic Neuromodulation, and Sensory Drive.

Following ref. 12, we modeled the time course of the mass cortical MUA using a simple dynamical systems model. Our model is a discrete-time generalization of the well-known FitzHugh–Nagumo (FHN) model (13). In brief, the model (schematized graphically in a simple form in Fig. 2*A* and with more detail in Fig. S4, and denoted hereafter as the “true-LC model”) describes cortical activity in an effective (rather than biophysically detailed) way. It contains terms of self-excitation and inhibition that lead to ongoing intrinsic fluctuations of excitability, an input from both i- and c-LC that models the potential causal effect of LC firing onto cortical activity, and an input term that describes the stimulus drive following the application of a somatosensory stimulus. This model, although simple, includes all of the three ingredients of which we want to describe the interactions (ongoing cortical state changes, neuromodulatory changes, and sensory drive).

The time course of mass cortical activity was modeled by the following finite-difference equations:**S1** and **S2**, and is effectively the convolution of

The function **S1** is defined as follows:

The term FHN is defined as follows:

Here, as in the spirit of ref. 12, the discrete time FHN model is considered just a phenomenological rather than biophysical description of cortical activity. We note, however, that this model lends itself to a mechanistic interpretation (12) as a model of mean field network dynamics. In fact

The term

The *SI Results* for our understanding of how these coupling terms affect the dynamics of our model). Although our current work did not determine whether this is the optimal way to include LC in cortical dynamics, we note that, as we shall report below, the nonlinear coupling term that depends on the *C*).

As it can be seen from Eq. **S5**, we also included in the model a parameter *SI Results*, *Computation of the Optimal Time Lag*, in which we report that for this dataset it was optimal to shift LC MUA back in the past by

To model stimulus-evoked responses, we also considered a “stimulus drive” current (which was called “kick function” in ref. 12), denoted as **S3**. The function *j* = 1,…,5) and

To evaluate the specific effect of the time course of LC firing on the dynamics of synchronized cortical states, we contrasted the performance of the true-LC model (in which the time evolution of cortical MUA depended both on past cortical activity and past LC activity) with that of a second model, named “no-LC” model and equivalent to the one previously introduced in ref. 12, in which evolution of cortical MUA depended only on its own past but not on LC firing. Namely, in the no-LC model, we did not include the LC input term into the model (or in other words, we set to zero the coupling terms **S5** of the LC input term). Comparison between the models with true-LC and no-LC inputs was used in this study to infer the contribution of LC to the dynamics of the cortical S1 MUA.

In other control analyses (*SI Results*, *Control Analysis on the Importance of the Terms Coupling LC Activity to the Model Cortical Dynamics*, *Control Analysis of the Differential Role of i-LC and c-LC*, and *Control Analysis of the Effect of Increasing the Number of Model Parameters*), we introduced additional modeling conditions that we called “phantom-LC” models, where the LC input (either both LCs or only one of the two LCs, and either in both fit and/or prediction period or in one of the two periods only, depending on the analysis) in the FHN equations was taken from the LC MUA in another randomly selected trial rather than in the current trial. The specific phantom-LC models and their utility will be explained when applied in *SI Results*.

It should be noted that, contrarily to the original FHN model, which was proposed to explain the dynamics of the action potential, which can take both positive and negative values, we are studying the dynamics of the mass firing rate, which is a positive definite variable. For this reason, we restricted the dynamical variable to be always positive during the whole dynamics. Moreover, because there is a physical limit for the maximum range of the firing rate of a population, we put an upper bound for

### Fitting the Model to the Data and Determining the Cortical State Parameters.

Here, we explain how we chose the model parameters that best fit the time course of the mass cortical MUA.

Because our discrete-time dynamical system (Eqs. **S1** and **S2**) is nonlinear and also depends on the neuromodulatory inputs from LC, the search for optimal parameters could be time-consuming because of the complicated input-dependent structure of the phase space, and therefore multiple possible local minima for the optimization problem could exist. Because of these reasons, we implemented a two-step procedure to fit the cortical activity to the model. In the first step, we used the recorded MUA to compute the right-hand side of Eq. **S1**, namely we defined **S7** and **S8** are linear with respect to the model parameters, and the optimization of the fit error between **S1** and **S2**). To avoid falling into local minima, we used a global optimization method, namely a genetic algorithm with large number of populations (a set of

To fit the model to the dynamics of a given stretch of single-trial cortical MUA, we minimized the error between the *dt* = 5-ms time resolution that was used to solve the dynamical system (Eqs. **S1** and **S2**). The window length for model fitting was chosen to be 1.5 s because dynamic cortical states are known to vary over timescales of seconds, and thus these windows are best suited to obtain accurate snapshots of ongoing dynamics (12).

The words “cortical states” and “cortical state variables” are used in the literature to describe a wide range of phenomena, with intrinsically different timescales. Following the work in ref. 12, we refer to “dynamic state” of cortex as to properties of network dynamics on a timescale of seconds or more (e.g., the synchronized and desynchronized states that can be detected from the local field potential spectrum). We use the term “activity state” to describe cortical states that persist on shorter timescales of few tens to few hundred milliseconds. Using this convention, the dynamical variables

### Model Integration and Quantification of Model Performance in the Prediction Period.

After having determined the model parameters by best fit to the cortical MUA during a 1.5-s window of spontaneous activity, we computed the cortical MUA that the model would predict after the end of the fit period. We generated a prediction of the cortical MUA dynamics in a postfit period by integrating numerically the model dynamics (Eqs. **S1** and **S2**) with the fit parameters using initial conditions given by the activity at the end of the fit period, for a further period (until 500 ms after the end-of-fit period when predicting continuation of spontaneous activity, and until 300 ms after the end of paw stimulation when predicting responses to stimulation), and still including the LC MUA contribution. Quantification of how well the dynamical model can predict cortical activity in the postfit period is useful because the test of model on data not used in the fitting period assures the cross-validation of the fitting quality.

For prediction of periods of spontaneous activity, we computed model performance simply as the Pearson correlation coefficient between the real MUA and its model prediction.

For prediction in periods of stimulus-induced activity, we computed two different quantifications of model performance.

The first measure quantifies the ability of the model to capture the variations in responses across stimuli and is computed as the correlation across different stimuli between the trial-averaged responses of model and data to each stimulus. We called this first measure (whose results are reported in Fig. S5*B*) the “signal correlation” between model and data, because signal correlation is a widely used term to refer to correlation between mean response profiles to stimuli (17, 28).

The second measure—which we call the noise correlation between the trial-to-trial variability of model and data—quantifies the ability of the model to predict the trial-to-trial variability of the real cortical MUA, and is defined as the correlation between the variation in each trial of real and model cortical MUA around their trial-averaged response to the considered stimulus. This quantity was computed separately for each stimulus and then averaged across all stimuli. We called this second measure (whose results are reported in Fig. 2*C* and Fig. S5 *G* and *H*) the “noise correlation” between model and data, because noise correlation is a widely used term to refer to correlation between variations of activity around the mean response to the considered stimulus (17, 28).

All these computations were performed for each time point of the postfit period. Then to quantify and present the results we took the mean correlation coefficients between 100 and 200 ms after the end of fitting period for the spontaneous activity trials (Fig. S5 *A*, *C*, *E*, and *F*) or after the end of the stimulus offset for the stimulation-evoked trials (Fig. 2*C* and Fig. S5 *G* and *H*). This period was chosen because, as we reported in Fig. 4*B*, this is the time period in which the information about the stimuli carried by cortical MUA reaches its peak, so this appeared a good choice for examining poststimulus responses.

To understand the extent to which the prediction of the single-trial variations around the stimulus-specific mean of the cortical stimulus-evoked response are stimulus independent and thus can be considered as an additive state-dependent variability term, we also computed for each trial the prediction of the variation around the stimulus mean of the single-trial cortical MUA (noise correlation) using the stimulus drive function from another stimulus different from the one that was actually applied in that trial. Then we computed the noise correlation between the trial-to-trial variability of model and data obtained when using the stimulus drive of another stimulus. Because the trials and stimulus drive functions used in this analysis were completely cross-validated and independent, there was no possible overfitting in these trial-to-trial noise correlation values. The result of this analysis is shown in Fig. S5*H* for the true-LC and no-LC models. Even though these noise correlation values were smaller than those obtained when using the stimulus drive function for the correct stimulus, they were still relatively large and statistically significant (*P* < 0.001, paired *t* test). Importantly, the noise correlation obtained with the stimulus drive function of another stimulus was still larger for the true-LC model than for the no-LC model, suggesting that the effect of the LC on the prediction of trial-to-trial variability is not affected by the fitting procedure of the stimulus drive function. The ability of the model to predict the trial-to-trial variability even with kick functions taken from other stimuli was crucial to the development of the algorithm that uses model prediction to discount state variability and increase the information that can be extracted from cortical MUA, because this observation allowed us to obtain single-trial model estimates of variability of cortical MUA that did not use any knowledge of the stimulus that was applied in each trial (*Discounting Trial-to-Trial Variability and Information Calculation*). An implication of the decrease in performance when using another stimulus drive function rather than the one corresponding to the stimulus actually presented in that trial suggests that the trial-to-trial variability originating from state dynamics and neuromodulation cannot always be described as purely additive, as proposed in some statistical models of state dependence (17).

### Time Domain Analysis of the Temporal Relationship Between LC and Cortical MUA.

To investigate whether the temporal order of firing of LC and MUA indicates the existence of a causal relationship between these variables, we performed an analysis of the relationship between the time course of i-LC and cortical MUA.

To quantify whether changes in i-LC firing tends to preferentially precede before or after, we computed the cross-correlogram (CCG) *i* in Eq. **S10** is performed over the *N* time points of the time series. With the sign convention used in Eq. **S10**, a CCG peak at positive (respectively negative) signifies that i,c-LC tends to fire before (respectively, after) cortex.

To measure the amount of directed transfer of information between the i-LC and cortical MUA, we computed the transfer entropy (TE) (10) between i,c-LC and cortex as the mutual information between the present values of cortical MUA and the past values of i,c-LC MUA, conditioned on the past values of cortical MUA. In other words, TE computes how well we can predict the present of cortical MUA when we know the past i,c-LC activity, after discarding the knowledge of current cortical MUA that can be obtained already from its own past. We note that TE attempts to measure causation in the Wiener–Granger sense, and it can be seen as a more general extension of Granger causality in that (being based on mutual information) it captures all possible kinds of relationships between the variables. In practice, however, conditioning on many past values of the neurophysiological signals makes it very difficult to sample the probabilities needed to compute mutual information. Following ref. 10, the empirical solution we chose is to use only one time delay, but to take it at a variable delay *B* for TE values corresponding to individual time delays), and then the delay value giving maximal delay for each electrode was chosen to plot results in Fig. 1*D*.

The TE between cortical MUA and i-LC MUA (quantifying whether cortex exerts a casual inference on i_LC firing in the Wiener–Granger sense) was computed in an analogous way as follows:

### Computation of Power Spectra of LC and Cortical MUA and of Their Correlation.

To investigate the relationship between the oscillatory structure of cortical and LC firing in the same trial, we computed single-trial power spectra of both cortical and LC MUA. The power spectrum of each trial was computed using the multitaper technique. The tapers are the discrete prolate spheroidal or Slepian, sequences defined in terms of their length *A*) was obtained as the difference between the power spectrum of the real MUA and the power spectrum of the model MUA normalized over the average of the power of real MUA between 1 and 10 Hz. For comparison, we plotted the power spectra of cortical and LC MUA without (Fig. 1*C*) and after (Fig. S6*A*) smoothing the real-data time series with the 100-ms half-Hanning window described above.

To compute the correlation between the power spectra of LC and cortical MUA, in each trial the LC MUA power was normalized by dividing it by the time average of the LC MUA in that trial. The values of correlation between LC power and real (respectively, true-LC model) cortical MUA power are reported in Fig. 1*F* (respectively, Fig. S7*C*).

### Computation of Synchronization Index.

To report concisely the range of cortical states observed in our dataset, we computed the synchronization index (S.I.), a widely used and very simple measure of brain state and anesthesia level (12, 17, 20). The S.I. was calculated as the ratio between the power of cortical MUA in the range of 1–5 Hz and the power in the range of 1–50 Hz. The distribution of S.I. over all of the experimental trials in all sessions is reported in Fig. S6*B*. This shows that most S.I. values were distributed in the range of 0.5–0.8, meaning that (as also apparent when visually examining the traces in Fig. 1 and Fig. S1) our experimental recordings were obtained from synchronized cortical networks.

### Phase Analysis.

To compute the phase of the cortical MUA at any instant of time, we first bandpassed the signal in the selected frequency band with a zero-lag Butterworth filter, we then computed the Hilbert transform, and we finally took its angle. For results reported in Fig. 1*E*, we took the 1- to 4-Hz delta cortical MUA band. For results reported in Fig. S3 *B* and *C*, we bandpassed the cortical MUA into 1-Hz-wide bands with center frequency from 1 to 10 Hz.

We computed the phase coherence between the model output and the real cortical MUA using the method explained in equations 8 and 9 of ref. 30.

To compute the phase dependence of the cortical and LC MUA, we first binned the cortical phase values into 16 equispaced divisions of the 2π cycle, and we averaged all MUA values (either of i-LC or of cortex) in each phase bin. We then computed and plotted the fraction of MUA elevation as function of phase, this fraction being defined as the value of MUA in a given phase bin minus the minimum value of MUA across all phase bins, divided by the average MUA. Thus, a value of for example 0.8 of this index would correspond to an elevation in firing rate with respect to the minimal phase whose value is 80% of the average MUA. These quantities were plotted in Fig. 1*E* and Fig. S3 *B* and *C*.

The results in Fig. 1*E* show the mean across the dataset and the bootstrap-estimated SD, whereas the results in Fig. S3 *B* and *C* plot the mean over the dataset.

### Computation of Correlation Between Fit Error and LC Power.

As we observed in the main text, the true-LC model described better the cortical MUA during the fit period (significantly smaller absolute fit error; Fig. 2*B*) than to the no-LC model. We asked whether the source of this fit quality improvement is related to the rhythmicity of LC MUA. To evaluate which oscillatory component of the LC MUA contributes more to the improvement of the true-LC model with respect to the no-LC model, we computed in each trial the difference between the absolute fit errors of the true-LC model and that of the no-LC model, and we then computed the correlation between this fit error difference and the i-LC power in different bands. If a particular band of i-LC MUA contributes more to helping the true-LC model to better explain the cortical MUA, we expect to observe negative correlation coefficient values between this fit error difference and the normalized power of the i-LC MUA in that band. We computed the power of LC for different 1-Hz-wide bands with center frequency between 1 and 10 Hz. Then we computed the correlation of the fit error differences of the two models with LC power. The results of this analysis are plotted in Fig. S7*B* and discussed in *SI Results*, *Additional Analyses of the Correlation of LC Firing Rate and LC Low-Frequency Power with Cortical Dynamics and Model Performance*.

### Discounting Trial-to-Trial Variability and Information Calculation.

As we reported above, the dynamical model predictions were able to predict—at least in part—the trial-to-trial variability in the cortical responses. We used the model prediction of the variation around the trial average of the neural response in a particular trial to “denoise” the neural response in that trial and thus increase the amount of information about stimuli that can be extracted by cortical MUA. The denoising of the single-trial cortical MUA was done by simply subtracting the prediction of the model of the expected variation around the mean in that trial. In the following, we describe the details of the denoising of neural responses and of the information calculation.

The first step of this calculation was to compute from the model the precise quantitative prediction of the variation around the mean of the MUA in a given trial. The dynamical systems model (true-LC or no-LC) gives for each instant of time a prediction of the stimulus-evoked cortical MUA in a single trial, given the choice of a stimulus drive function corresponding to one of the stimuli in the set. It is important to note that computation of this prediction using the stimulus drive function of the stimulus actually presented in that trial would require knowing which stimulus was presented. However, using stimulus knowledge for computing stimulus information may increase information artificially, introducing sources of information other than the stimulus information content of neural responses. We thus decided to generate—for the information analysis—an additive prediction of single-trial variation of the cortical MUA around the mean that does not need the knowledge of the stimulus. Our work reported above (Fig. S5*H*) showed that the model generates the best prediction of variability when using the stimulus drive function with the same stimulus, but it still gives a good prediction even when using the stimulus drive function of another stimulus. Therefore, to generate a denoising that does not use knowledge of the stimulus (and thus can be implemented by a real online decoder of neural responses), we computed the model prediction using a single stimulus drive function (corresponding to one of the stimuli in the experimental set) for all of the trials. We chose the stimulus drive function of the stimulus that maximized the so-obtained information.

We then refined these model predictions of variations of responses around the stimulus mean that do not rely on the knowledge of the stimulus that we actually presented by rescaling them into the final prediction of the amount of MUA to be subtracted in each trial using a cross-validated linear regression. We first obtained the cross-validated linear regression coefficient between the actual and the predicted variations of responses around the mean removing that trial (“training set”), and then we rescaled the output of the model in the considered trial with this linear regression coefficient. This was our estimate of the model prediction of the MUA to be subtracted in each trial. We note that we found that the rescaled cross validation had little practical effect on results, but we found it useful to avoid any possible spurious variability reduction or information increase due to either overfitting in the variability discounting procedure.

We also estimated how well one could denoise cortical MUA based on the following simple proxies of the state of cortex at the time of stimulus presentation: the S.I. during the state detection period, the delta-band phase of cortical MUA (indicated as

We computed for each such proxy of cortical state the prediction of the variation around the mean to be subtracted in each trial using the same cross-validated linear regression as for the LC model estimates as detailed above. In short, we first obtained the cross-validated linear regression coefficient between the actual variations of responses around the mean and the proxy state parameter in the same trial by first removing the considered trial (training set), and then we multiplied the value of the state proxy in the current trial by the this linear regression coefficient to obtain our estimate of the MUA to be subtracted in that trial.

We then computed the information that these denoised cortical MUA values carry about the somatosensory stimuli, as explained below.

It is important to bear in mind that the above prediction of state-dependent variability was obtained by the model in a cross-validated way and without using any knowledge about which stimulus was presented, and that the subtraction of estimated variability may not always be the optimal way to discount it. It is possible that future studies will reveal more effective ways to reduce state variability by using procedures more refined than a simple subtraction. Thus, the information values that we obtain in this study after we discount by subtraction the estimate of the trial-to-trial variability should be taken as robust lower bounds to the information that can be gained by discounting the state-dependent variability captured by the models.

### Computation of Information from MUA Before and After Discounting the Trial-to-Trial Variability Predicted by the Model.

For a given quantification

By evaluating the information carried by different ways of defining and quantifying the response

The computation of information requires the estimation of the stimulus–conditional response probabilities. These probabilities are not known a priori and must be measured experimentally from a finite number of trials. The estimated probabilities suffer from finite sampling errors, which induce a systematic error (bias) in estimates of the information (28). In this article, unless otherwise stated, we used the following procedure to correct for the bias when computing information from real data. First, to facilitate the sampling of its probability, we considered responses *r* that were binned into three bins. Although the discretization facilitates the sampling of neural response probabilities, the information measures still suffer from limited sampling bias. We thus used a Panzeri–Treves (28) to estimate and subtract the bias of each information quantity.

To check that the so-obtained information estimates were not biased, we computed a “permuted” information estimates from data after shuffling at random the association between stimuli and responses, and we always found that the so-obtained information were two orders of magnitude smaller than the information in real data during the response windows with maximal information (observed in the 125- to 200-ms post–end-of-stimulus window; Fig. 4*B*). This check on permuted data confirmed that our information estimates do not suffer from an upward bias.

All of the information computations were done using the information toolbox in Matlab (28).

### Simulation of the Effect of Phasic LC Stimulation on Cortical Responses to Stimuli.

We used our model to study how phasic activation of LC may shape the cortical information representation of sensory stimuli. Phasic activation of LC refers to brief transient bursts of firing of LC neurons, which may happen spontaneously or be brought about by electrical stimulation.

We simulated the effect of phasic LC activation as follows. We created a true-LC model with the best-fit parameters from each of the trials in our dataset, and then we applied to each of these models a stimulus drive function generating responses to the different stimuli that we again took from fits to real data. For each of these trials, we then run to different set of computer simulations. In the first set, describing normal spontaneously occurring states, we run the model using the true-LC time course in that trial as neuromodulatory input. In the second set, we simulated the effect of phasic activation of LC by adding to the true-LC dynamics a transient moment of firing whose amplitude and timing were chosen to maximize the amplitude resulting response (Fig. 5*A*). We then inferred the specific effect of phasic LC activation by comparing across these two simulations the cortical stimulus-evoked MUA.

The burst of phasic LC activity was considered to have duration of 20 ms (modeling a phasic activation) with amplitude (A) and time (T) before the stimulus onset of the burst being parameters that were optimized as follows. We examined a range of amplitudes, from zero to 5 SDUs of the LC MUA, and a range of T from 0 to 200 ms before the stimulus onset. The response of each phasic burst of stimulation was considered as the mean MUA in the late 125- to 200-ms post–end-of-stimulus component of the stimulus-evoked response (the response period presenting maximal sensory information in real data) after the end of the stimulus offset. For each stimulus, we considered the optimal pair of parameters (A_{opt},T_{opt}) as the phasic LC pulse that produces the maximum response.

For the plot in Fig. 5*B*, we sorted the stimulus rank according to the response values according to their mean MUA stimulus-evoked activity in the real data in the 125- to 200-ms post–end-of stimulus period computed across all recordings.

In the main text (Fig. 5*C*), we reported the values of peak time and rise time of the stimulus-evoked response predicted by the models with and without phasic LC activation. In these results, the peak time was computed as the time when the response gets its maximum value in the range 0–300 ms after end of stimulus. The rise time was computed as the earliest time before the peak time in which the cortical MUA is increasing and reaches to 50% of the peak response.

## SI Results

### Computation of the Optimal Time Lag.

To obtain the optimal model parameter *D*) that there was a clear modulation of the fit error with LC time lag and an optimal time lag equal to 20 ms (implying that the LC had the maximal effect on cortical activity with a 20-ms delay). We used this optimal value of time lag in all reported analyses of the true-LC model.

The difference between the fit error at zero lag and the fit error with a lag of 20 ms was significant (*P* < 0.01, paired *t* test). The result that the true-LC model fits cortical MUA better when shifting the LC activity in the past rather than considering LC activity with no delay or shifted toward the future suggests that in these data the temporal modulation of LC firing exerts a causal effect on the state of cortical firing, rather than LC firing being driven by cortical state dynamics.

### Values of the Best-Fit Model Parameters over the Population.

To gain further insights into how the coupling between LC and cortex may operate in our model, we quantified the values of the best-fit parameters of the cortical–LC coupling observed across experiments. Results are reported in Fig. S8.

We discuss here in particular the coupling parameters

### Investigation into How the Nonlinear Coupling Between LC and Cortex May Increase the Power of Model Cortical Low-Frequency Oscillations.

In the following, we present our intuitive reasoning about why the LC input, and in particular its nonlinear coupling term, can generate additional low-frequency oscillation power in the model cortical MUA. For simplicity, we focus the discussion only on the effect of the i-LC, which is the one contributing the most to the cortical state (Fig. S5 *E* and *G*), and we neglect the effect of c-LC. As it can be seen from **Eqs. S1–S3**, the total LC input **S5**) that depends both on the product of LC MUA and of cortical MUA. For obtaining a total LC input that can change sign over time, it is necessary that the coefficients of the linear and nonlinear terms have opposite signs. Although we did not impose any constraint on the signs of the input coefficients when we fitted the model to the data, we indeed found (Fig. S8 *G* and *I*) that the best-fit parameters had opposite signs for *E*)] would give rise to a less negative nonlinear term than if the same LC burst happened at the peak of cortical excitability. The net effect observed in data and model is that overall LC bursts increase the power of low-frequency cortical fluctuations (Fig. 3*A*). To understand why the net effect of LC was to increase the cortical power, and to characterize better what could be the effect of LC at different phases of the up-down cortical cycle, we investigated how the total LC current of the model depended on the cortical delta-band (1–4 Hz) phase. Given that we decide to focus only on the effect of the i-LC, we first computed the LC current **S5**. To subtract out the oscillatory component that the ups and downs of cortical MUA add to this current and to focus specifically on the pure contribution of i-LC firing to the input current, we subtracted from

The dependence of the pure-LC current on the phase of cortical MUA in the 1- to 4-Hz delta-band is shown in Fig. S9*A* both for the true-LC model (constructed taking as input LC firing from the same trial; red line) and for the phantom-LC model (constructed taking as input LC firing from a random trial; blue line). We found that the pure-LC current for true-LC model had a positive sign in the ascending phase up to just before the cortical MUA peak, and it was negative in the descending phase after the peak of cortical excitability. In comparison, the pure-LC current for the phantom-LC model showed no dependence on phase (Fig. S9*A*), suggesting that the timing of LC burst is crucial to produce the asymmetry between ascending and descending phase of the up-down cycle in the sign of the LC net effect.

To further understand how the pure-LC current can produce oscillations coherent with cortical MUA and can increase its power, we considered a linear first-order approximation of the contribution of the pure-LC current input to the dynamics of cortical MUA by solving the following decoupled part of Eq. **S1**:*B* the result of computing **S14**, as function of the phase of the delta-band cortical MUA. The results of this analysis show that, on average over all of the real data, the specific phase coupling between the LC bursts and cortical MUA produces an oscillatory component that is coherent in phase with the cortical MUA and hence can increase its power. All in all, this analysis suggests that properly timed LC bursts are crucial to increase the power of cortical low-frequency fluctuations because they produce a different and phase-coherent effect on the ascending and descending phase of cortical up- and down-state changes.

### Control Analysis on the Importance of the Terms Coupling LC Activity to the Model Cortical Dynamics.

To better understand the importance of the nonlinear coupling between LC and cortex, we examined the effect of modifying the LC input term defined in Eq. **S5**. We tested the role of the

We found that having the *P* < 0.05, paired *t* test). Importantly, we also found that having the *C* show that the true-LC model (that contains both the *P* < 0.05, paired *t* test) better than the model without

### Control Analysis of the Differential Role of i-LC and c-LC.

To understand which of the two LCs contributed the most to cortical dynamics, we created models that selectively took out the contribution of one of the two LCs in the same trial, and we compared the performance of models that did or did not consider the contribution of each LC. To remove the contribution of one of the two LCs, we first replaced its LC time series with the MUA taken from the same LC in another randomly selected phantom trial and then we performed the model fitting. The additional models that we so created were denoted in this paper and the axis labels of Figs. S5 and S7 as follows: the (*i-T*,*c-P*) model is the one that takes as input the true i-LC firing from the same trial and the phantom c-LC firing from a random trial; the (*i-P*,*c-T*) model is the one taking as input the phantom i-LC firing from a random trial and the true c-LC firing from the same trial; the (*i-P*,*c-P*) model is the one taking as input the firing from a random trial for both LCs. The comparison of these two models is useful to reveal which LC may contribute more to cortical S1 MUA dynamics.

The results of the prediction quality for these models are reported in Fig. S5*E* (prediction of spontaneous activity) and in Fig. S5*G* (prediction of stimulus-evoked response). In both cases, the (*i-T*,*c-P*) model containing only the contribution of the i-LC performed nearly as well as the true-LC model containing the contribution of both LCs, and it performed much better than the (*i-P*,*c-T*) model containing the contribution of the c-LC. The power fit error of different models (Fig. S7*A*) showed a similar trend.

Taken together, these results show that, according to our model, the i-LC had a much stronger impact on the dynamics of cortical MUA than the one exerted by the c-LC, consistent with anatomical observations (9). These results are supported by the model-free TE analysis reported in *Supplemental Results on the Time Domain Analysis of the Temporal Relationship Between LC and Cortical MUA*.

### Control Analysis of the Effect of Increasing the Number of Model Parameters.

We also used the (*i-P*,*c-P*) model that takes as input the firing from a random trial for both LCs to investigate whether the increase in performance of the true-LC model over the no-LC model was simply due to the true-LC model having more parameters, or whether it rather reflected the fact that the insertion of true-LC activity in the same trial genuinely increase the ability of the model to describe cortical MUA. We found (Fig. S5 *E* and *G*) that across the dataset the (*i-P*,*c-P*) model performed always significantly (*P* < 0.01, paired *t* test) less well than the true-LC model. The performance of the (*i-P*,*c-P*) model was comparable (*P* > 0.1, paired *t* test) to that of the no-LC model. These results show that the increase in performance of the true-LC model reflected the ability of this model to capture a real influence of LC on the dynamics of cortical MUA, and it is not only a trivial consequence of the increase in the number of model parameters.

### Control Analyses on the Differential Role of LC Firing in Improving State Detection and in Improving Stimulus Drive.

LC firing can influence the accuracy of the model’s postfit prediction in two ways: either by permitting to obtain more accurate “dynamic-state variables” in the fit (state detection) period, or by contributing a more accurate instantaneous drive during the postfit period. To evaluate which of these contributions is more important, we performed the following control analyses. When we replaced only in the fit period the simultaneous LC activity taken from the considered trial with LC firing taken from a phantom random trial (thereby destroying only the contribution of LC to state detection), we obtained a large and significant (*P* < 0.01, paired *t* test) reduction of the model’s prediction of trial-to-trial response variability (Fig. S5*F*). When we replaced LC firing with a randomly selected phantom period of LC activity only in the prediction period (thereby destroying only the LC drive during the prediction period but leaving state detection intact), we obtained a smaller and not significant (*P* > 0.5, paired *t* test) reduction of the model’s prediction of trial-to-trial response variability (Fig. S5*F*). These results suggest that a major contribution of LC firing to the prediction of cortical dynamics comes from help determining the correct dynamic-state variables.

### Additional Analyses of the Correlation of LC Firing Rate and LC Low-Frequency Power with Cortical Dynamics and Model Performance.

To better understand the relationship between LC firing rhythmicity and cortical dynamics, we performed the following analyses. In the main text, we reported that the power of low-frequency LC rhythmic activity correlates with the power of low-frequency cortical MUA. To test how specific this effect was with regard to the rhythmicity of LC MUA, we investigated the correlation between the time-averaged firing rate of LC and the power of cortical MUA oscillations. We found that the delta cortical power did not correlate (*P* > 0.2, paired *t* test) with LC averaged firing rate (Fig. S3*A*), and thus the amplification of cortical state variations was specifically related to low-frequency rhythmicity of LC activity.

To better understand how critical was the low-frequency rhythmicity of LC activity to the ability of the model to reproduce better cortical dynamics, we performed a further analysis (Fig. S7*B*), in which we found that the model’s error fit was anticorrelated to the trial-to-trial variations in the power of the LC low-frequency rhythmic activity, but was not correlated to the power of the LC activity at higher frequencies, suggesting again that the low-frequency LC rhythmicity has a pivotal role in predicting the low-frequency cortical dynamics.

### Comparing the Prediction of Trial-to-Trial Response Variability Obtained with the Dynamical Systems Modeling to That Obtained with Simpler State Indices.

To further understand whether the models predict the trial-to-trial variability of cortical responses to stimuli better than other simple quantifications of the network state or excitability around the time of stimulus application, we computed the noise correlation between the variation of the poststimulus MUA around its mean and the following three “proxies” of cortical state: the MUA at the time point just preceding the stimulus application, the delta-phase (1–4 Hz) of the cortical MUA at the moment of stimulus onset, and the value of the S.I. in the prestimulus period (obtained as the ratio of the power in the 1–5 Hz and the power in the 1–10 Hz). Reassuringly, we found that both dynamical system models of cortical state outperformed by 85–180% these reference performance values obtained with the simple state proxies (Fig. 2*C*). The better predictability of trial-to-trial noise correlation of the true-LC and no-LC dynamical system models with respect to the three simple proxies of cortical state translated also into a larger information advantage when using these predictions of trial-to-trial variability to discount state variations from cortical MUA and gain sensory information (Fig. S10*A*).

### Control Analyses for the Information Calculation.

We performed several additional analyses to both better understand the importance of specific aspects of LC dynamics into the performance of the state variability discount to increase the sensory information that can be extracted from cortical MUA, or to check that these discounts could not be due to artifacts.

We first asked what could be the differential role of c-LC and i-LC in achieving a better state variability discount and a better information extraction from cortical MUA. To investigate this question, we computed the information after discounting the state variability using (as described in *Control Analysis of the Differential Role of i-LC and c-LC*) models in which one of the two LCs in the considered trial with another period of activity from the same LC in another randomly selected period. We found (Fig. S10*B*) that, on average over the 125- to 200-ms post–end-of-stimulus period, randomly replacing the c-LC did not lead to a significant (*P* > 0.5, paired *t* test) decrease of the model’s ability to gain information by discounting state variability, whereas randomly replacing the i-LC led to a significant (*P* < 0.05, paired *t* test) decrease. These results further confirm that the dominant contribution to cortical state changes comes from the i-LC.

To evaluate if LC MUA influences the stimulus information gain by contributing to obtain more accurate “state parameters” in the fit period, or by the information that it might itself carry about the stimulus, we performed another control analysis. Given two sets of stimuli, 1 and 2, at each poststimulus time *i* (*B* and show that the informationless-LC model has the same performance in extracting information as the true-LC model. This means that the higher information in the true-LC model was entirely due to its better prediction of cortical state, and rules out that the increased information of the true-LC model was due to information contained in the poststimulus LC firing rate.

### Control Analyses for the Effect of Jitter in the Time of LC Phasic Activation on the Trial-to-Trial Reliability of the Temporal Profile of Cortical Responses to Stimuli.

The results presented in the main text (*Results*, *Modeling the Effect of LC Phasic Activation on Cortical Responses to Stimuli*) show that a phasic activation of the LC that precedes in a finely timed way the application of a sensory stimulus increases the trial-to-trial reliability of the timing of neural responses to the sensory stimuli. How accurately timed does the LC activation preceding the stimulus need to be to reduce response variability?

To address this question, we repeated the simulations of the addition of a LC phasic activation onto the spontaneous LC firing by choosing an addition time that was randomly jittered around the optimal time value by a jitter value distributed in the interval [−X/2 +X/2] ms, where X quantified the extent of the allowed time jitter and was varied parametrically in range of 10–200 ms across different simulations. (We always used the optimal phasic LC amplitude in all of the time jitter values.) We found (Fig. S11) that we obtained a significant (*P* < 0.001, paired *t* test) reduction of trial-to-trial in response peak time and rise time during the late response window when using phasic LC activation as long at the jitter range was below 120 ms for peak time and below 40 ms for rise time.

### Supplemental Results on the Time Domain Analysis of the Temporal Relationship Between LC and Cortical MUA.

We computed the CCG between i-LC MUA and cortical MUA, and we found that the results were compatible with the view that i-LC firing precedes changes of cortical MUA. On average over the entire dataset, CCG was higher (*P* < 0.001, paired *t* test) when i-LC MUA was shifted back in time in the range of up to 25 ms (Fig. S2*A*). The peak of the CCG when LC activity was shifted back in time was significantly larger than the value of CCG at time 0 for three out of four animals at *P* < 0.001 (paired *t* test).

Moreover, we computed TE [a well-known nonlinear information-theoretic statistical measure of causality in the Wiener–Granger sense (10)] between the time series of i-LC and cortical MUA. As shown in Fig. 1*D*, on average over the entire dataset TE from i-LC to S1 cortex was approximately three times larger than TE in the opposite direction with *P* < 0.001 (paired *t* test). The TE from i-LC to cortex was significantly larger than the TE in the opposite direction for three out of four animals at *P* < 0.001 (paired *t* test).

The TE from c-LC to S1 cortex (Fig. 1*D* and Fig. S2*B*) followed a similar pattern to the TE from i-LC to S1 cortex, but it was smaller in magnitude (*P* < 0.01, paired *t* test). This suggests, that, in accordance with anatomical considerations (9) and other analyses reported in this study (*SI Results*, *Control Analysis of the Differential Role of i-LC and c-LC*), i-LC is the LC influencing ipsilateral cortex the most.

These results are all compatible with the hypothesis that LC firing rapidly exerts a causal influence on cortex.

## Acknowledgments

We thank J. Assad, T. Fellin, P. Series, and members of S.P.’s laboratory for useful feedback. This work was supported by the European Commission (FP7-ICT-2011.9.11/284553, “SICODE,” and FP7-2007-2013/PITN-GA-2011-290011, “ABC”), by the Max Planck Society, and by the Autonomous Province of Trento (“Grandi Progetti 2012,” “ATTEND”).

## Footnotes

- ↵
^{1}To whom correspondence may be addressed. Email: stefano.panzeri{at}iit.it or nikos.logothetis{at}tuebingen.mpg.de.

Author contributions: H.S., O.E., N.K.L., and S.P. designed research; H.S., R.N., O.E., and S.P. performed research; H.S. and S.P. contributed new reagents/analytic tools; H.S. and S.P. analyzed data; and H.S., R.N., O.E., N.K.L., and S.P. wrote the paper.

The authors declare no conflict of interest.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1516539112/-/DCSupplemental.

Freely available online through the PNAS open access option.

## References

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵.
- Eschenko O,
- Magri C,
- Panzeri S,
- Sara SJ

- ↵
- ↵
- ↵
- ↵.
- Curto C,
- Sakata S,
- Marguet S,
- Itskov V,
- Harris KD

- ↵.
- Fitzhugh R

*Bull Math Biophys*17:257–278 - ↵
- ↵.
- Schölvinck ML,
- Saleem AB,
- Benucci A,
- Harris KD,
- Carandini M

- ↵
- ↵
- ↵
- ↵
- ↵.
- Kayser C,
- Wilson C,
- Safaai H,
- Sakata S,
- Panzeri S

- ↵
- ↵
- ↵
- ↵
- ↵.
- Lemieux M,
- Chen JY,
- Lonjers P,
- Bazhenov M,
- Timofeev I

- ↵.
- Usher M,
- Cohen JD,
- Servan-Schreiber D,
- Rajkowski J,
- Aston-Jones G

- ↵.
- Braintenberg V,
- Schuetz A

- ↵
- ↵.
- Marzo A,
- Totah NK,
- Neves RM,
- Logothetis NK,
- Eschenko O

- ↵

## Citation Manager Formats

## Article Classifications

- Biological Sciences
- Neuroscience

## Sign up for Article Alerts

## Jump to section

## You May Also be Interested in

### More Articles of This Classification

### Biological Sciences

### Neuroscience

### Related Content

### Cited by...

- Resting brain fluctuations are intrinsically coupled to visual response dynamics
- The locus coeruleus is a complex and differentiated neuromodulatory system
- Catecholaminergic Manipulation Alters Dynamic Network Topology Across Behavioral States
- Phasic locus coeruleus activity regulates cortical encoding of salience information
- Surprise About Sensory Event Timing Drives Cortical Transients in the Beta Frequency Band
- Back to Pupillometry: How Cortical Network State Fluctuations Tracked by Pupil Dynamics Could Explain Neural Signal Variability in Human Cognitive Neuroscience
- Multiple Transient Signals in Human Visual Cortex Associated with an Elementary Decision
- Catecholaminergic Neuromodulation Shapes Intrinsic MRI Functional Connectivity in the Human Brain