Skip to main content
  • Submit
  • About
    • Editorial Board
    • PNAS Staff
    • FAQ
    • Accessibility Statement
    • Rights and Permissions
    • Site Map
  • Contact
  • Journal Club
  • Subscribe
    • Subscription Rates
    • Subscriptions FAQ
    • Open Access
    • Recommend PNAS to Your Librarian
  • Log in
  • My Cart

Main menu

  • Home
  • Articles
    • Current
    • Special Feature Articles - Most Recent
    • Special Features
    • Colloquia
    • Collected Articles
    • PNAS Classics
    • List of Issues
  • Front Matter
  • News
    • For the Press
    • This Week In PNAS
    • PNAS in the News
  • Podcasts
  • Authors
    • Information for Authors
    • Editorial and Journal Policies
    • Submission Procedures
    • Fees and Licenses
  • Submit
  • About
    • Editorial Board
    • PNAS Staff
    • FAQ
    • Accessibility Statement
    • Rights and Permissions
    • Site Map
  • Contact
  • Journal Club
  • Subscribe
    • Subscription Rates
    • Subscriptions FAQ
    • Open Access
    • Recommend PNAS to Your Librarian

User menu

  • Log in
  • My Cart

Search

  • Advanced search
Home
Home

Advanced Search

  • Home
  • Articles
    • Current
    • Special Feature Articles - Most Recent
    • Special Features
    • Colloquia
    • Collected Articles
    • PNAS Classics
    • List of Issues
  • Front Matter
  • News
    • For the Press
    • This Week In PNAS
    • PNAS in the News
  • Podcasts
  • Authors
    • Information for Authors
    • Editorial and Journal Policies
    • Submission Procedures
    • Fees and Licenses

New Research In

Physical Sciences

Featured Portals

  • Physics
  • Chemistry
  • Sustainability Science

Articles by Topic

  • Applied Mathematics
  • Applied Physical Sciences
  • Astronomy
  • Computer Sciences
  • Earth, Atmospheric, and Planetary Sciences
  • Engineering
  • Environmental Sciences
  • Mathematics
  • Statistics

Social Sciences

Featured Portals

  • Anthropology
  • Sustainability Science

Articles by Topic

  • Economic Sciences
  • Environmental Sciences
  • Political Sciences
  • Psychological and Cognitive Sciences
  • Social Sciences

Biological Sciences

Featured Portals

  • Sustainability Science

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
Research Article

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

View ORCID ProfileHouman Safaai, Ricardo Neves, Oxana Eschenko, Nikos K. Logothetis, and Stefano Panzeri
PNAS October 13, 2015 112 (41) 12834-12839; first published September 28, 2015; https://doi.org/10.1073/pnas.1516539112
Houman Safaai
aNeural Computation Laboratory, Istituto Italiano di Tecnologia, 38068 Rovereto, Italy;
  • Find this author on Google Scholar
  • Find this author on PubMed
  • Search for this author on this site
  • ORCID record for Houman Safaai
Ricardo Neves
bDepartment of Physiology of Cognitive Processes, Max Planck Institute for Biological Cybernetics, 72076 Tübingen, Germany;
  • Find this author on Google Scholar
  • Find this author on PubMed
  • Search for this author on this site
Oxana Eschenko
bDepartment of Physiology of Cognitive Processes, Max Planck Institute for Biological Cybernetics, 72076 Tübingen, Germany;
  • Find this author on Google Scholar
  • Find this author on PubMed
  • Search for this author on this site
Nikos K. Logothetis
bDepartment of Physiology of Cognitive Processes, Max Planck Institute for Biological Cybernetics, 72076 Tübingen, Germany;
cDivision of Imaging Science and Biomedical Engineering, University of Manchester, Manchester M13 9PT, United Kingdom
  • Find this author on Google Scholar
  • Find this author on PubMed
  • Search for this author on this site
  • For correspondence: stefano.panzeri@iit.it nikos.logothetis@tuebingen.mpg.de
Stefano Panzeri
aNeural Computation Laboratory, Istituto Italiano di Tecnologia, 38068 Rovereto, Italy;
  • Find this author on Google Scholar
  • Find this author on PubMed
  • Search for this author on this site
  • For correspondence: stefano.panzeri@iit.it nikos.logothetis@tuebingen.mpg.de
  1. Contributed by Nikos K. Logothetis, August 20, 2015 (sent for review June 22, 2015)

  • Article
  • Figures & SI
  • Info & Metrics
  • PDF
Loading

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.

  • oscillations
  • state dependence
  • dynamical systems
  • information coding
  • somatosensation

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. 1C) 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. 1C) 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. 1B 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. 1B and Fig. S1 C and D).

Fig. 1.
  • Download figure
  • Open in new tab
  • Download powerpoint
Fig. 1.

Relationship between firing of LC and cortex. (A) Example trial showing the time course of data and model during the state detection and the prediction period of spontaneous activity. (Top) Raster of spikes recorded in each cortical channel. (Middle) Mass cortical MUA (sum of spikes from all cortical channel, plotted as black dashed line and expressed in units of spikes per second per channel) is compared with the output of true-LC (red line) and no-LC (green) model. (Lower) Time course of i-LC (black line) and c-LC (gray line) MUA, plotted in SD units (SDUs). (B) Example trial showing data and model during the state detection period and the prediction period of stimulus-evoked activity. Conventions are as in A, with the addition of a series of gray vertical pulses plotting the model’s stimulus-drive function obtained in this trial. (C) Power spectra (normalized by the total power in the 1- to 50-Hz range) of cortical MUA (Left) and i-LC MUA (Right). (D) TE computed from either c- or i-LC to S1 or in the opposite direction. Power spectra (normalized by the total power in the 1- to 50-Hz range) of cortical MUA (Left) and i-LC MUA (Right). Results plotted as mean ± SEM over the dataset. (E) The fraction of elevation (with respect to the phase in which MUA is minimal) of the real MUA of i-LC and cortex plotted against the phase of cortical activity in delta (1–4 Hz)-band. Results are plotted as mean ± SEM across dataset (n = 830 spontaneous trials) in C and as mean across dataset and bootstrap-estimated SD in D and E. (F) Average over the dataset of the Pearson correlation between i-LC power (normalized by the i-LC firing rate in the trial) and cortical power at each frequency.

Fig. S1.
  • Download figure
  • Open in new tab
  • Download powerpoint
Fig. S1.

More example trials. (A and B) Two example trials showing the time course of real and of model MUA during the fit period of spontaneous activity (Left) and during a continued prediction period of spontaneous activity (Right). Conventions are as in Fig. 1A. (C and D) Two example trials showing the time course of real and of model MUA during the fit period of spontaneous activity (Left) preceding a somatosensory stimulus and during the prediction period that includes activity recorded after the application of the somatosensory stimulus (Right). Conventions are as in Fig. 1B.

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. S2A). 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. 1D; see also Fig. S2). Together, these results support the hypothesis that LC firing exerts a causal influence on cortex within few tens of milliseconds.

Fig. S2.
  • Download figure
  • Open in new tab
  • Download powerpoint
Fig. S2.

Cross-correlogram (CCG) and transfer entropy (TE). (A) The result shows the CCG between the i,c-LC and S1 MUA time series as a function of the time lag Δ entering Eq. S10 (with the convention that positive time lags indicate LC preceding cortical S1 activity). (B) TE from cortex to i-LC (black) and c-LC (gray) and from i-LC and c-LC to cortex is plotted as function of the time lag Δ entering its definition (Eqs. S11 and S12). In this figure, results are plotted as mean over the dataset ± bootstrap-estimated SD (n = 830 spontaneous trials).

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. 1F) 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. S3A), 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).

Fig. S3.
  • Download figure
  • Open in new tab
  • Download powerpoint
Fig. S3.

Additional results on temporal relationships between firing of LC and of cortex. (A) The correlation coefficient between the time average of the i-LC firing rate over a 1.5-s spontaneous activity trial and the cortical MUA power in the same window (mean ± bootstrap-estimated SD over dataset). (B) Color plot of the fraction of elevation (with respect to the phase in which MUA is minimal) of real i-LC MUA as function of the phase (y axis) of the bandpassed cortical MUA frequency (x axis) during spontaneous activity. (C) Color plot of the fraction of elevation (with respect to the phase in which MUA is minimal) of real cortical MUA as function of the phase (y axis) of the bandpassed cortical MUA frequency (x axis) during spontaneous activity. B and C report results averaged over the dataset.

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. 1E). 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. 2A 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).

Fig. 2.
  • Download figure
  • Open in new tab
  • Download powerpoint
Fig. 2.

Schematics and performance of the model. (A, Top) Time course of cortical MUA is modeled as a self-exciting and inhibiting dynamical system coupled to a neuromodulatory input (MUA recorded in the LCs) and a stimulus drive. (Bottom) The model was fitted to real cortical MUA during a state detection period of spontaneous activity (from t = −1.5 to 0 s), and then it was solved in the prediction period (from t = 0 s onward). (B) Time average during the state detection period of the normalized absolute fit error of the no-LC and true-LC model. (C) Time average over the 100- to 200-ms range of the prediction period of the noise correlation between the variation around the trial average of cortical MUA at a given time and the following five predictions: MUA predicted by either no-LC or true-LC model, by synchronization index (S.I.) during the state detection period, by delta-band phase of cortical MUA (indicated as ϕ0) or by the MUA amplitude (ν0) at stimulus onset. In B and C, graphs plot mean ± SEM over the dataset (n = 830 spontaneous trials, n = 826 evoked trials). In all figures, the single asterisk (*), double asterisk (**), and triple asterisk (***) indicate significance values of P < 0.05, 0.01, and 0.001, respectively.

Fig. S4.
  • Download figure
  • Open in new tab
  • Download powerpoint
Fig. S4.

Details of dynamical model. The figure schematizes the details of different components of the discrete-time version of the FHN dynamical model that we used. The external input currents to the FHN model (the LC input containing the contribution of i- and c-LC, the stimulus drive current, and the tonic current) are shown on the Left, whereas the FHN variables are plotted on the Right and the coupling coefficients are plotted in the Middle.

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. 1B 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. 2B) 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. S5A). 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. S5B; 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. 2C; 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. S6C) or the cortical delta-phase or MUA amplitude at the time of stimulus application (Fig. 2C), demonstrating the value of dynamical systems in describing cortical states.

Fig. S5.
  • Download figure
  • Open in new tab
  • Download powerpoint
Fig. S5.

Control analyses to determine the role of specific model components. (A) Time average over the 100- to 200-ms range of the prediction period of correlation between the cortical MUA and its model predictions for spontaneous activity. (B) Time average over the 100- to 200-ms range of the prediction period of the signal correlation across different stimuli of the mean response to each stimulus of real data and models. (C) Time average over the 100- to 200-ms range of the prediction period of the correlation between the cortical MUA and different model predictions for the spontaneous activity. The models compared are as follows: the true-LC model containing both the linear (ci−LC,cc−LC) and nonlinear (di−LC,dc−LC) LC coupling terms; the “true-LC (c)” model containing only the linear (ci−LC,cc−LC) LC coupling terms; the no-LC model not containing any coupling to LC. (D) The average over time during the state detection period of the normalized absolute fit error for model predictions with different values of the Δ parameter quantifying the lag between LC firing and its effect over cortical MUA. (E) Time average over the 100- to 200-ms range of the prediction period of the correlation between the cortical MUA and different model predictions for the spontaneous activity. The models compared here are as follows: the true-LC model, which is the model that takes as input the true i-LC and c-LC firing from the same trial; the (i-T,c-P) model is the model 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 model that takes 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 model that takes as input the firing from a random trial for both LCs; the no-LC model is the model that does not take any LC input. (F) Time average over the 100- to 200-ms range of the prediction period of the correlation between the cortical MUA and different model predictions for the spontaneous activity. The models compared here are as follows: the (pr-T,po-T) is the model that takes in the model the true i-LC and true c-Lc firing from the same trial both for the state detection period (referred here as “pr” for pre) and the prediction period (referred here as “po” for post); the (pr-T,po-P) is the model that takes in the model the true i-LC and true c-Lc firing from the same trial for the state detection period but takes the phantom i-LC and phantom c-Lc firing from a random trial for the prediction period; the (pr-P,po-T) is the model that takes in the model the phantom i-LC and phantom c-Lc firing from a random trial for the state detection period but takes for the prediction period the true i-LC and true c-Lc firing from the same trial; the (pr-P,po-P) is the model that takes in the model the phantom i-LC and phantom c-Lc firing from a random trial both for the state detection and prediction periods. (G) Time average over the 100- to 200-ms range of the prediction period of the noise correlation between the variation around the trial average of the stimulus-evoked cortical MUA at a given time and its prediction using different models with different combinations of true and phantom LC inputs (see E for nomenclature of models). (H) Time average over the 100- to 200-ms range of the prediction period of the noise correlation between the variation around the trial average of the stimulus-evoked cortical MUA at a given time and its prediction using different models using either the same kick-function for each stimulus or a different kick-function from the original stimulus of each trial. In this figure, all bars graphs show mean ± SEM over the dataset (n = 830 spontaneous trials, n = 826 evoked trials). In all figures, the single asterisk (*), double asterisk (**), and triple asterisk (***) indicate significance values of P < 0.05, 0.01, and 0.001, respectively.

Fig. S6.
  • Download figure
  • Open in new tab
  • Download powerpoint
Fig. S6.

Power spectrum and distribution of the values of the synchronization index. (A) Plot (mean ± SEM over the dataset of n = 830 spontaneous trials) of the normalized power of the smoothed cortical MUA (Left) and smoothed i-LC MUA. The normalization was done over the total power in the range of 1–50 Hz in each trial. The cortical MUA signal and LC signals were smoothed with a causal Henning window of 100 ms. (B) The fraction of spontaneous trials (n = 830) showing a given synchronization index (S.I.). Results show the distribution over all trials of 1.5 s of spontaneous activity that we used to fit the model to the data.

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. 3A), 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. S7C) 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).

Fig. 3.
  • Download figure
  • Open in new tab
  • Download powerpoint
Fig. 3.

Model performance in the frequency domain. (A) Difference between power spectrum of real MUA and MUA of true-LC and no-LC models, normalized by the mean power of real MUA power in the 1- to 10-Hz band. The black horizontal line shows the frequency range with significant difference between the two models. (B) The phase coherence between the delta-bands of the real cortical MUA and of the MUA predicted by each of the two models. A and B plot mean ± SEM over the dataset (n = 830 spontaneous trials).

Fig. S7.
  • Download figure
  • Open in new tab
  • Download powerpoint
Fig. S7.

Model ability of predicting temporal structure of LC–cortex relationship. (A) The mean ± SEM over the dataset (n = 830 spontaneous trials) of the difference between the power spectrum of the 1- to 10-Hz real cortical MUA and the power spectrum in the same range of the output of each different models using different combinations of true and phantom LC as input (see Supporting Information and Fig. S5E for model’s nomenclature). Results are normalized over the mean real-MUA power between 1 and 10 Hz. (B) The correlation coefficient between the difference in fit error between the true-LC and the no-LC model and the i-LC power (curves plot mean ± bootstrap-estimated SD over dataset, n = 830 spontaneous trials). (C) Average over the dataset (n = 830 spontaneous trials) of the Pearson correlation between the i-LC power (normalized to the i-LC firing rate) in a given trial at a given frequency and power of the true-LC model’s cortical MUA at each frequency.

Fig. S8.
  • Download figure
  • Open in new tab
  • Download powerpoint
Fig. S8.

Model parameters for different LC input conditions. (A–J) The mean ± SEM over the dataset (n = 830 spontaneous trials, n = 826 evoked trials) of the value of each of the model’s parameters for the true-LC (red) model that includes the contribution of both LCs and the phantom-LC (blue) model that is fit to the data after taking as input the time series of the LC MUA taken from another randomly selected trial. Each panel shows the value of a different model parameter. The parameters (except τ) are computed from the model of Eqs. S1 and S2 using the cortical MUA computed in the units of spikes per millisecond and the LC MUA in SD units, exactly as explained in Supporting Information. The τ parameter reported in F is in milliseconds after multiplying the unitless τ obtained from the finite-difference model Eqs. S1 and S2 with dt=5 ms.

Fig. S9.
  • Download figure
  • Open in new tab
  • Download powerpoint
Fig. S9.

Role of the LC current in increasing the model cortical power. (A) The pure-LC current (see Supporting Information for definition) quantifying the pure effect of LC on cortex is computed either using the true-LC model (red line) and the phantom-LC (blue line), and is plotted against the cortical MUA phase in delta-band. (B) The MUA obtained from integrating the pure-LC current (Eq. S14) is plotted against the cortical MUA phase in delta-band for true-LC (red line) and phantom-LC (blue line). The results in this figure show the mean ± SEM over the dataset (n = 830, spontaneous trials).

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. 3B), 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. 4A, 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. 4A) 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.

Fig. 4.
  • Download figure
  • Open in new tab
  • Download powerpoint
Fig. 4.

Gaining sensory information by discounting model predictions of trial-to-trial variability. (A) This mechanism is exemplified with a single trial taken from real data. Stimulus 1 was presented in this trial, and the cortical MUA response in this trial (dark-blue full line) is plotted against the trial-averaged response to stimulus 1 (dashed blue) and to stimulus 2 (dashed red) for a 0- to 300-ms period after stimulus offsets. The true-LC model predicts that the response variation around the trial average in this trial was positive (thus, the model prediction of what should be subtracted—shown with black arrows—was negative). The subtraction of the predicted variability from the single-trial response produces a “discounted” response (light-blue full line) that is much closer to the trial-averaged response to the stimulus 1 that was presented in that trial, and it is thus much easier to decode, than the original response. (B) Time course of mutual information about the somatosensory stimuli computed from either the original cortical MUA at a given time, or after discounting from it the trial-to-trial variability using the predictions of each model. Results are plotted as mean ± SEM over n = 43 stimulus pairs.

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. 4B. 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. S10A). This shows the worth of dynamical systems for extracting information from neural data.

Fig. S10.
  • Download figure
  • Open in new tab
  • Download powerpoint
Fig. S10.

Information computation controls. (A) Mutual information about the somatosensory stimuli carried either by the original cortical MUA at a given time (black bar), or computed after discounting from the original cortical MUA in each trial the estimate of its variation around the trial average obtained using the following models: the true-LC model, the no-LC model, the MUA ν0 at the time point just preceding the stimulus application, the delta-phase ϕ0 of delta-band (1–4 Hz) cortical MUA at the moment of stimulus application, the value of the S.I. in the prestimulus period. (B) Mutual information about the somatosensory computed after discounting in each trial from the original cortical MUA the estimate of its variation around the mean obtained using the following models: the true-LC model that takes as input the true i-LC and c-LC firing from the same trial; the (i-T,c-P) model 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 that takes 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 that takes as input the firing from a random trial for both i-LC and c-LC; the no-LC model, which does not take any LC input; and a further “informationless-LC” control model in which we added to the poststimulus activity a noise term canceling out any information that the LC firing may carry about the stimulus. In both panels, results are reported as mean ± SEM over the dataset (n = 43 stimulus pairs).

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. 5A) 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. 5B) 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. 5C). 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).

Fig. 5.
  • Download figure
  • Open in new tab
  • Download powerpoint
Fig. 5.

Simulating the effect of phasic LC input on the model stimulus responses. (A) Schematic of the simulated experiment: A phasic LC activation (shown in purple) is added to the recorded spontaneous LC firing (shown in gray). The stimulus-evoked MUA of the model after phasic stimulation is then compared with the stimulus-evoked MUA of the same model without the addition of the phasic LC pulse. (B) The cortical MUA model response to each of the stimuli averaged in the 125- to 200-ms post–end-of-stimulus response period. Full lines report stimulus-evoked responses with (purple) and without (gray) LC phasic activation, whereas the dashed black line reports spontaneous activity averaged in a randomly selected 75-ms period before the stimulus. For this plot, the nine stimuli indicated on the x axis were sorted according to their mean response across all recordings. (C) The SD across trials at fixed stimulus of the model’s peak response time and of the rise time, with and without addition of the phasic LC activation to the model. Throughout this figure, results are plotted as mean ± SEM over the dataset (n = 826 evoked trials).

Fig. S11.
  • Download figure
  • Open in new tab
  • Download powerpoint
Fig. S11.

The effect of phasic LC input timing on the model stimulus responses properties. (A) The SD across trials at fixed stimulus of the model’s peak response time with (purple line) and without (black line) addition of a time jitter (whose range is plotted on the x axis) between the time of simulated LC phasic activation and of simulated application of the somatosensory stimulus. (B) Similar to A but for the model’s response rise time. In both panels, time jitters were applied around the optimal delay between simulated phasic LC activation and application of simulated somatosensory stimulus that was used in Fig. 5. Throughout this figure, results are plotted as mean ± SEM over the dataset (n = 826 evoked trials).

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 νMUA, computed as the convolution of the spike train with half-Hanning causal window of 100-ms width (i.e., the value of νMUA in the ith 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 ν in the LC ipsilateral (i-LC) and contralateral (c-LC) to the recorded S1 cortex was computed as the absolute value of the bandpassed (400–3,000 Hz, fourth-order Butterworth) time series of the extracellular potential. This time series was then smoothed exactly as explained above for the cortical MUA. For the purpose of plotting it in Fig. 1 and Fig. S1, the MUA from the LCs was converted into SD units (SDUs), although this normalization does not affect the results of the model fitting. The MUAs from the LCs were used as an input to some of the models of cortical MUA (described below) and were indicated in the model equations (see below) as νi−LC and νc−LC, respectively.

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. 2A 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:νmodel(t+dt)=νmodel(t)+f(νmodel(t),wmodel(t),Iext(t)),[S1]wmodel(t+dt)=wmodel(t)+(νmodel(t)−wmodel(t))τ.[S2]In the above, time derivatives are computed by one–time-step differentiation, and νmodel and wmodel are the phase-space variables of the dynamics of the discrete-time FHN model. The model variable νmodel represents the mass cortical MUA (whose experimental counterpart is the experimental cortical MUA, νMUA described above). The slower model variable wmodel is defined by the finite-difference equations Eqs. S1 and S2, and is effectively the convolution of νmodel with an exponential kernel of time constant τ (τ is one of the model parameters determined by best fit). It describes the integrated past activity of the dynamical system on a timescale given by time constant τ. We use a time step of dt=5 m to numerically solve the model’s equations.

The function f in the model equation Eq. S1 is defined as follows:f(νmodel(t),wmodel(t),Iext(t))=FHN+ILC(t)+Iconst(t)+ks(t),[S3]where different terms are defined and discussed in the following.

The term FHN is defined as follows:FHN=∑i=13ai(νmodel(t))i+bwmodel(t).[S4]FHN is a standard polynomial term used in FHN model (13). The FHN model is a well-studied simplification of the Hodgkin–Huxley model of action potential generation in single neurons. However, the dynamical equations of the model are also used to study a wide range of oscillatory dynamical systems, including the neural mass models describing the mesoscopic-scale interaction of neural systems in ref. 12. This is because the FHN model presents both fast (ν) and slow (w) variables that can interact in an inhibitory and excitatory way (including self-excitation and self-inhibition) and that (as we shall see below) have an interpretation in terms of effective description of neural population variables.

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 νmodel(t) can be seen as the total rate of excitatory and inhibitory cells of the network and wmodel(t) as the combined effects of adaptive phenomena such as synaptic depression and cellular-level accommodation that reduce network excitability. The model parameters ai represent recurrent excitation and inhibition, and b relates to the extent to which increases in wmodel reduce the excitability of the network (12).

The term Iconst is a constant tonic drive on all cells, which may, e.g., reflect a tonic increase in mean input firing rate from the thalamus, or any other source of tonic and unspecific firing.

The ILC term is the term representing the total effect of firing of i-LC and c-LC on cortical activity and is defined as follows:ILC(t)=(ci−LC+di−LCνmodel(t))νi−LC(t−Δ)+(cc−LC+dc−LCνmodel(t))νc−LC(t−Δ).[S5]This term expresses the neuromodulatory inputs as a sum of a first “linear” coupling term (the one with coupling coefficients ci−LC and cc−LC) that is purely linear in terms of LC activity, and a second “nonlinear” coupling term (the one with coupling coefficients di−LC and dc−LC) that depends on the product between the MUA of LC and of cortex. This kind of additional linear coupling term expands the dynamical range of the possible solutions (for example, negative values of di−LC or dc−LC can generate slower modes in the dynamics; see 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 di−LC and dc−LC coefficients significantly improves the agreement between model and cortical data (Fig. S5C).

As it can be seen from Eq. S5, we also included in the model a parameter Δ representing the relative delay with which LC firing affects cortical dynamics. The value of Δ is a model parameter that should be determined empirically. We determined the optimal Δ delay value by minimizing the normalized absolute fit error over the dataset for the true-LC model (see 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 Δ=20 ms).

To model stimulus-evoked responses, we also considered a “stimulus drive” current (which was called “kick function” in ref. 12), denoted as ks in Eq. S3. The function ks is zero during periods of spontaneous activity and is nonzero at the time when the sensory stimuli are applied. In the case of the dataset that we considered, each stimulus consisted of a set of five pulses with a particular frequency and amplitude. We modeled the stimulus drive function of a given stimulus s as a sequence of alpha functions each placed at the time of each of single pulses. Thus, each kick function corresponding to a particular stimulus was written as follows:ks(t)=∑j=15αj(t−tj)et−tjβ,[S6]where tj is time of the onset of the j th pulse, and αj (j = 1,…,5) and β are model parameters. The six parameters of the stimulus drive function were optimized—separately for each of the nine different somatosensory stimuli used in the experiments—by minimizing the fit error (computed over all of the trials in a particular stimulus in a given session) between the cortical MUA of model and of real data model in the period extending from the stimulus onset until 300 ms after the end of the stimulation. (To avoid contamination of MUA with stimulation artifacts, the extracellular signals recorded from both the cortical and LC channels were set to zero for 3 ms each side of each stimulation pulse.)

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 ci−LC, cc−LC, di−LC, and dc−LC in the Eq. 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 νmodel(t) equal to the twice of maximum value of νMUA during the whole experiment and the firing was discharged whenever it reaches to the upper bound. This was to prevent divergent solutions during the optimization and also to prevent the inputs (LC and stimulus inputs) to push the solution out of its plausible dynamical range.

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 νself from solving the following finite-difference equations:νself(t+dt)=νself(t)+f(νMUA(t),wMUA(t),Iext(t)),[S7]wself(t+dt)=wself(t)+(νMUA(t)−wself(t))τ.[S8]The right-hand sides of Eqs. S7 and S8 are linear with respect to the model parameters, and the optimization of the fit error between νself and νMUA can be done with a simple linear regression. Note that, in this case, the neuromodulatory input also does not depend on the model dynamical parameter νself because LC activity is coupled with νMUA. In the second step, we used the solutions that we obtained for the “self-dynamics” as the initial values for the full dynamical equations (Eqs. 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 106 chromosome populations). This two-step optimization procedure was empirically found to increase considerably the quality and speed of the fitting procedure.

To fit the model to the dynamics of a given stretch of single-trial cortical MUA, we minimized the error between the νMUA(t) measured from cortex and its model prediction νmodel(t), the error being defined for each trial as the normalized absolute-value error:NAE=∑t|(νmodel(t)−νMUA(t))|∑tνMUA(t),[S9]where in the above expression the sums are over the time points spanning the full fitting period of length 1.5 s during spontaneous activity. Time points were summed with the same time interval 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(t) and wmodel(t)—as well as νi−LC(t) and νc−LC(t)—describe the activity state of the network at any instant, whereas the model parameters Θ=[ai,b,ci−LC,di−LC,cc−LC,dc−LC,τ] that specify how νmodel(t) and wmodel(t) evolve in time describe its dynamic state.

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. S5B) 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. 2C 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. 2C and Fig. S5 G and H). This period was chosen because, as we reported in Fig. 4B, 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. S5H 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) R(Δ) between the time series of i,c-LC MUA νi,c−LC(t) and of cortical MUA νMUA(t) as a function of the time lag Δ between the two time series, as follows:R(Δ)=1N∑t=1N(νMUA(t)−ν¯MUA(t))(νi,c−LC(t−Δ)−ν¯i,c−LC(t−Δ))σMUAσi,c−LC,[S10]where ν¯i,c−LC and σi,c−LC are the mean and SDs of the i-LC, and ν¯MUA and σMUA are the mean and SDs of the cortical MUA. The sum over 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 Δ, which is the same for both time series and whose value is varied parametrically within a range to test the potential effect of causations at different delays. As a result, TE from i,c-LC to cortical MUA was computed as follows:TEνi,c−LC→νMUA=I(νMUA(t);νi,c−LC(t−Δ)|νMUA(t−Δ)),[S11]where in the above I denotes mutual information between the two variables. To control the biases entering into the probability estimates needed in computing the entropies, we first binned the signals into 5 equispaced bins (we checked that 3–10 bins give similar results, but we used the one which produced more robust results in our dataset), and then we computed the corresponding entropies using the information toolbox (28). As stated above, the delay Δ between past and present was varied parametrically (in the range of 5–75 ms; see Fig. S2B 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. 1D.

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:TEνMUA→νi,c−LC=I(νi,c−LC(t);νMUA(t−Δ)|νi,c−LC(t−Δ)).[S12]

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 L and their bandwidth W. For each choice of these parameters, there are K=2LW−1 tapers highly concentrated in frequency and can be averaged for spectral estimation. The optimal choice of the parameters should be done with empirical considerations. We took LW=2, which appeared to give good trade-off between the resolution of the power spectrum estimation and providing sharper dependencies of the power to the frequency. We checked the results were robust for changing the value of LW between 1 and 4. The normalized power spectrum fit error (reported in Fig. 3A) 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. 1C) and after (Fig. S6A) 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. 1F (respectively, Fig. S7C).

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. S6B. 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. 1E, 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. 1E and Fig. S3 B and C.

The results in Fig. 1E 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. 2B) 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. S7B 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. S5H) 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 ϕ0), and the MUA amplitude (indicated as ν0) at the time of stimulus onset. Note that the S.I. used in this specific analysis was calculated as the ratio of the power in the range of 1–5 Hz over the power in the range of 1–10 Hz because we found empirically that this range was the one that maximized the correlation with the trial-to-trial variability. Other ranges for definitions of S.I. (for example, using the 1- to 50-Hz range for the denominator), yielded slightly less strong correlations with the trial-to-trial variability.

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 r of a neural response, we computed the information carried by about the stimulus can be quantified by Shannon’s Mutual Information equation (14), abbreviated hereafter as Information:I(S;R)=∑s,rP(s)P(r|s)log2P(r|s)P(r),[S13]where P(s) is the probability of presentation of stimuluss, P(r|s) is the probability of observing response r given presentation of stimulus s, and P(r) is the probability of observing response r across all trials to any stimulus. This equation quantifies the maximal reduction of uncertainty (i.e., the gain in information) about the stimuli obtained from a single-trial observation of a neuronal response (averaged over all stimuli and responses). Information is measured in bits (1 bit corresponds to a reduction of uncertainty by a factor of 2) and is an upper bound on the amount of knowledge about stimuli that can be extracted by any decoding algorithm operating on neural responses r.

By evaluating the information carried by different ways of defining and quantifying the response r, we can evaluate the capacity of different candidate neural codes. Here, we computed two types of information. The first was the information carried by the poststimulus MUA, which corresponds to the traditional calculation of information performed without knowledge about the state that led to this response. The second was the information carried by the MUA after subtracting out the estimate of trial-to-trial variability obtained with various predictions methods that were described above.

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. 4B). 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. 5A). 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 (Aopt,Topt) as the phasic LC pulse that produces the maximum response.

For the plot in Fig. 5B, 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. 5C), 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 Δ of the true-LC model quantifying the time lag after which the LC firing affects the cortical MUA, we computed the normalized absolute mean fit error for models that had true-LC inputs of different time lags with respect to cortical MUA. The time lags that we used to search for optimal parameters ranged between −100 and +100 ms. We found (Fig. S5D) 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 ci,c−LC and di,c−LC between the cortical MUA and the LCs ipsilateral and contralateral to the cortex S1 we recorded from. Notably, we found that the ipsilateral di−LC parameter was negative. This finding was specific to the ipsilateral coupling (the dc−LC contralateral coupling term was not different from zero) and it was genuine, because it became zero when coupling the cortical model to “phantom” LC firing taken from another randomly selected trial (phantom-LC model results in Fig. S8). As we will discuss in the next section, the fact that di−LC was negative and ci−LC was positive is important to understand how the model may account for compatible with the empirical observations reported in the main text that suggest that LC firing contributes to increase the amplitude of slow cortical oscillations.

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 ILC is one of the terms contributing to the rate at which the MUA changes at each time step. A positive (respectively, negative) total LC input term ILC at time t tends to increase (respectively, decrease) the cortical MUA at time t+dt. The reason why the linear LC input term is not enough—by itself—to create additional oscillatory components to the MUA can be understood intuitively as follows. If the LC total input ILC contained only the linear term defined as ci−LCνi−LC with no other nonlinear term, given that MUA is nonnegative the sign of the change over time of cortical MUA would be determined by the sign of the coefficient ci−LC. Thus, its contribution would not change sign over time and would be unable to create additional oscillatory modes to the cortical MUA. To create additional oscillatory components, it is necessary to have an additional nonlinear input term that can change the sign of the total LC input term ILC during time. The simplest way to achieve a time-changing total LC input is to add to the linear component of the LC coupling a second nonlinear component di−LCνi−LCνMUA (Eq. 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 ci−LC and di−LC (positive sign for ci−LC and negative sign for di−LC). It is apparent that, if ci−LC>0 and di−LC<0, the LC input term ci−LCνi−LC+di−LCvi−LCνMUA can be both positive and negative depending on the relative values of the νi−LC and νMUA at each time. In particular, with this coupling term the contribution of a phasic LC burst to the model’s dynamics would depend on the cortical state at which it happens. If the LC burst happens during the down state of cortical MUA, which implies that the absolute value of nonlinear coupling term is small, the effect of the LC to the dynamics would be described by a purely additive term that pushes cortex toward a more excitable state. If the LC burst happens at a time of elevated cortical firing, then the total LC input can be positive or negative depending on the relative strength of the cortical and LC MUA. For example, an LC burst happening a few tens of milliseconds before the peak of cortical excitability [as it happens more often than not (Fig. 1E)] 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. 3A). 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 Ii−LC from the i-LC MUA using the first term of Eq. 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 Ii−LC the “shifted current” I˜i−LC computed from the model best-fit parameters after shifting the i-LC MUA ν˜i−LC in each trial with a random time delay. The so-obtained current, defined as Ii−LC−I˜i−LC, includes only the part of the LC current that is modulated by the LC rhythmic bursts and by their phase coupling with the cortical MUA, and was called “pure-LC current” hereafter.

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. S9A 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. S9A), 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:Δνpure−LC−current(t+dt)−Δνpure−LC−current(t)=Ii−LC(t)−I˜i−LC(t).[S14]We report in Fig. S9B the result of computing Δνpure−LC−current both for true-LC (red line) and phantom-LC (blue line), after summing over Eq. 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 di,c−LC terms by analyzing models with and without these terms (in other words, in these latter models we fixed the nonlinear coupling parameters di,c−LC to have zero value).

We found that having the di,c−LC terms reduced significantly the fit error (P < 0.05, paired t test). Importantly, we also found that having the di,c−LC terms led also to better predictions of the cortical activity. Results presented in Fig. S5C show that the true-LC model (that contains both the ci,c−LC and di,c−LC coupling terms) predicts the trial-to-trial variability of spontaneous cortical MUA better (P < 0.05, paired t test) better than the model without di,c−LC coupling does, implying that the coupling terms improve the ability of the model to predict cortical dynamics. Note that the model without di,c−LC coupling does not predict trial-to-trail variability better than the no-LC model does, implying that the linear ci,c−LC coupling terms is not enough to model the effect of LC on cortex.

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. S5E (prediction of spontaneous activity) and in Fig. S5G (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. S7A) 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. S5F). 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. S5F). 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. S3A), 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. S7B), 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. 2C). 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. S10A).

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. S10B) 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 t we first computed the mean LC activity of each of the stimulus sets and then we construct an “informationless LC” by adding a noise term to the original LC signal that perfectly equalizes the LC trial-averaged LC MUA between the two stimuli. For each of the stimuli, this noise was sampled from a uniform random distribution with amplitude εi=2×(νLC¯−νLCi¯), where νLCi¯ is the mean of LC activity of stimulus i (i=1,2) and νLC¯ is the grand mean of the LC MUA across all of the trials to both stimuli. The so-obtained informationless-LC signal carries negligible information about the stimulus at each poststimulus time. Thus, when it is given as input to the model, it cannot artificially enhance the estimate of stimulus information from cortical MUA based on discounting the state variability from its model prediction. (Note also that, because this noise did not have any temporal structure, the resulting informationless-LC signal will have similar low-frequency structure as of the original LC signal.) We thus used the prediction of this informationless-LC model to predict the trial-to-trial variability and to discount it from the single-trial cortical MUA. Results of this control analyses are reported in Fig. S10B 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. S2A). 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. 1D, 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. 1D and Fig. S2B) 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

  • ↵1To 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.

View Abstract

References

  1. ↵
    1. Buonomano DV,
    2. Maass W
    (2009) State-dependent computations: Spatiotemporal processing in cortical networks. Nat Rev Neurosci 10(2):113–125
    .
    OpenUrlCrossRefPubMed
  2. ↵
    1. Douglas RJ,
    2. Martin KA
    (2007) Mapping the matrix: The ways of neocortex. Neuron 56(2):226–238
    .
    OpenUrlCrossRefPubMed
  3. ↵
    1. Harris KD,
    2. Thiele A
    (2011) Cortical state and attention. Nat Rev Neurosci 12(9):509–523
    .
    OpenUrlCrossRefPubMed
  4. ↵
    1. Lee SH,
    2. Dan Y
    (2012) Neuromodulation of brain states. Neuron 76(1):209–222
    .
    OpenUrlCrossRefPubMed
  5. ↵
    1. Sara SJ,
    2. Bouret S
    (2012) Orienting and reorienting: The locus coeruleus mediates cognition through arousal. Neuron 76(1):130–141
    .
    OpenUrlCrossRefPubMed
  6. ↵
    1. Berridge CW,
    2. Waterhouse BD
    (2003) The locus coeruleus-noradrenergic system: Modulation of behavioral state and state-dependent cognitive processes. Brain Res Brain Res Rev 42(1):33–84
    .
    OpenUrlCrossRefPubMed
  7. ↵
    1. Aston-Jones G,
    2. Cohen JD
    (2005) An integrative theory of locus coeruleus-norepinephrine function: Adaptive gain and optimal performance. Annu Rev Neurosci 28:403–450
    .
    OpenUrlCrossRefPubMed
  8. ↵
    1. Eschenko O,
    2. Magri C,
    3. Panzeri S,
    4. Sara SJ
    (2012) Noradrenergic neurons of the locus coeruleus are phase locked to cortical up-down states during sleep. Cereb Cortex 22(2):426–435
    .
    OpenUrlAbstract/FREE Full Text
  9. ↵
    1. Lee SB,
    2. Beak SK,
    3. Park SH,
    4. Waterhouse BD,
    5. Lee HS
    (2009) Collateral projection from the locus coeruleus to whisker-related sensory and motor brain regions of the rat. J Comp Neurol 514(4):387–402
    .
    OpenUrlCrossRefPubMed
  10. ↵
    1. Schreiber T
    (2000) Measuring information transfer. Phys Rev Lett 85(2):461–464
    .
    OpenUrlCrossRefPubMed
  11. ↵
    1. Aston-Jones G,
    2. Bloom FE
    (1981) Activity of norepinephrine-containing locus coeruleus neurons in behaving rats anticipates fluctuations in the sleep-waking cycle. J Neurosci 1(8):876–886
    .
    OpenUrlAbstract
  12. ↵
    1. Curto C,
    2. Sakata S,
    3. Marguet S,
    4. Itskov V,
    5. Harris KD
    (2009) A simple model of cortical dynamics explains variability and state dependence of sensory responses in urethane-anesthetized auditory cortex. J Neurosci 29(34):10600–10612
    .
    OpenUrlAbstract/FREE Full Text
  13. ↵
    1. Fitzhugh R
    (1955) Mathematical models of threshold phenomena in the nerve membrane. Bull Math Biophys 17:257–278
    .
  14. ↵
    1. Quian Quiroga R,
    2. Panzeri S
    (2009) Extracting information from neuronal populations: Information theory and decoding approaches. Nat Rev Neurosci 10(3):173–185
    .
    OpenUrlCrossRefPubMed
  15. ↵
    1. Schölvinck ML,
    2. Saleem AB,
    3. Benucci A,
    4. Harris KD,
    5. Carandini M
    (2015) Cortical state determines global variability and correlations in visual cortex. J Neurosci 35(1):170–178
    .
    OpenUrlAbstract/FREE Full Text
  16. ↵
    1. Goris RLT,
    2. Movshon JA,
    3. Simoncelli EP
    (2014) Partitioning neuronal variability. Nat Neurosci 17(6):858–865
    .
    OpenUrlCrossRefPubMed
  17. ↵
    1. Ecker AS, et al.
    (2014) State dependence of noise correlations in macaque primary visual cortex. Neuron 82(1):235–248
    .
    OpenUrlCrossRefPubMed
  18. ↵
    1. Sachidhanandam S,
    2. Sreenivasan V,
    3. Kyriakatos A,
    4. Kremer Y,
    5. Petersen CCH
    (2013) Membrane potential correlates of sensory perception in mouse barrel cortex. Nat Neurosci 16(11):1671–1677
    .
    OpenUrlCrossRefPubMed
  19. ↵
    1. Lecas JC
    (2004) Locus coeruleus activation shortens synaptic drive while decreasing spike latency and jitter in sensorimotor cortex. Implications for neuronal integration. Eur J Neurosci 19(9):2519–2530
    .
    OpenUrlCrossRefPubMed
  20. ↵
    1. Kayser C,
    2. Wilson C,
    3. Safaai H,
    4. Sakata S,
    5. Panzeri S
    (2015) Rhythmic auditory cortex activity at multiple timescales shapes stimulus-response gain and background firing. J Neurosci 35(20):7750–7762
    .
    OpenUrlAbstract/FREE Full Text
  21. ↵
    1. Vanrullen R,
    2. Busch NA,
    3. Drewes J,
    4. Dubois J
    (2011) Ongoing EEG phase as a trial-by-trial predictor of perceptual and attentional variability. Front Psychol 2:60
    .
    OpenUrlCrossRefPubMed
  22. ↵
    1. Thut G,
    2. Miniussi C,
    3. Gross J
    (2012) The functional importance of rhythmic activity in the brain. Curr Biol 22(16):R658–R663
    .
    OpenUrlCrossRefPubMed
  23. ↵
    1. Schroeder CE,
    2. Lakatos P
    (2009) Low-frequency neuronal oscillations as instruments of sensory selection. Trends Neurosci 32(1):9–18
    .
    OpenUrlCrossRefPubMed
  24. ↵
    1. Destexhe A
    (2011) Intracellular and computational evidence for a dominant role of internal network activity in cortical computations. Curr Opin Neurobiol 21(5):717–725
    .
    OpenUrlCrossRefPubMed
  25. ↵
    1. Lemieux M,
    2. Chen JY,
    3. Lonjers P,
    4. Bazhenov M,
    5. Timofeev I
    (2014) The impact of cortical deafferentation on the neocortical slow oscillation. J Neurosci 34(16):5689–5703
    .
    OpenUrlAbstract/FREE Full Text
  26. ↵
    1. Usher M,
    2. Cohen JD,
    3. Servan-Schreiber D,
    4. Rajkowski J,
    5. Aston-Jones G
    (1999) The role of locus coeruleus in the regulation of cognitive performance. Science 283(5401):549–554
    .
    OpenUrlAbstract/FREE Full Text
  27. ↵
    1. Braintenberg V,
    2. Schuetz A
    (1998) Cortex: Statistics and Geometry of Neuronal Connectivity (Springer, Berlin)
    .
  28. ↵
    1. Magri C,
    2. Whittingstall K,
    3. Singh V,
    4. Logothetis NK,
    5. Panzeri S
    (2009) A toolbox for the fast information analysis of multiple-site LFP, EEG and spike train recordings. BMC Neurosci 10:81
    .
    OpenUrlCrossRefPubMed
  29. ↵
    1. Marzo A,
    2. Totah NK,
    3. Neves RM,
    4. Logothetis NK,
    5. Eschenko O
    (2014) Unilateral electrical stimulation of rat locus coeruleus elicits bilateral response of norepinephrine neurons and sustained activation of medial prefrontal cortex. J Neurophysiol 111(12):2570–2588
    .
    OpenUrlAbstract/FREE Full Text
  30. ↵
    1. Mormann F,
    2. Lehnertz K,
    3. David P,
    4. Elger CE
    (2000) Mean phase coherence as a measure for phase synchronization and its application to the EEG of epilepsy patients. Physica D 144(3-4):358–369
    .
    OpenUrlCrossRef
PreviousNext
Back to top
Article Alerts
Email Article

Thank you for your interest in spreading the word on PNAS.

NOTE: We only request your email address so that the person you are recommending the page to knows that you wanted them to see it, and that it is not junk mail. We do not capture any email address.

Enter multiple addresses on separate lines or separate them with commas.
Modeling the effect of locus coeruleus firing on cortical state dynamics and single-trial sensory processing
(Your Name) has sent you a message from PNAS
(Your Name) thought you would like to see the PNAS web site.
CAPTCHA
This question is for testing whether or not you are a human visitor and to prevent automated spam submissions.
Citation Tools
Locus coeruleus firing and cortical state dynamics
Houman Safaai, Ricardo Neves, Oxana Eschenko, Nikos K. Logothetis, Stefano Panzeri
Proceedings of the National Academy of Sciences Oct 2015, 112 (41) 12834-12839; DOI: 10.1073/pnas.1516539112

Citation Manager Formats

  • BibTeX
  • Bookends
  • EasyBib
  • EndNote (tagged)
  • EndNote 8 (xml)
  • Medlars
  • Mendeley
  • Papers
  • RefWorks Tagged
  • Ref Manager
  • RIS
  • Zotero
Request Permissions
Share
Locus coeruleus firing and cortical state dynamics
Houman Safaai, Ricardo Neves, Oxana Eschenko, Nikos K. Logothetis, Stefano Panzeri
Proceedings of the National Academy of Sciences Oct 2015, 112 (41) 12834-12839; DOI: 10.1073/pnas.1516539112
Digg logo Reddit logo Twitter logo Facebook logo Google logo Mendeley logo
  • Tweet Widget
  • Facebook Like
  • Mendeley logo Mendeley
Proceedings of the National Academy of Sciences: 112 (41)
Table of Contents

Submit

Sign up for Article Alerts

Article Classifications

  • Biological Sciences
  • Neuroscience

Jump to section

  • Article
    • Abstract
    • Results
    • Discussion
    • Materials and Methods
    • SI Materials and Methods
    • SI Results
    • Acknowledgments
    • Footnotes
    • References
  • Figures & SI
  • Info & Metrics
  • PDF

You May Also be Interested in

Abstract depiction of a guitar and musical note
Science & Culture: At the nexus of music and medicine, some see disease treatments
Although the evidence is still limited, a growing body of research suggests music may have beneficial effects for diseases such as Parkinson’s.
Image credit: Shutterstock/agsandrew.
Large piece of gold
News Feature: Tracing gold's cosmic origins
Astronomers thought they’d finally figured out where gold and other heavy elements in the universe came from. In light of recent results, they’re not so sure.
Image credit: Science Source/Tom McHugh.
Dancers in red dresses
Journal Club: Friends appear to share patterns of brain activity
Researchers are still trying to understand what causes this strong correlation between neural and social networks.
Image credit: Shutterstock/Yeongsik Im.
White and blue bird
Hazards of ozone pollution to birds
Amanda Rodewald, Ivan Rudik, and Catherine Kling talk about the hazards of ozone pollution to birds.
Listen
Past PodcastsSubscribe
Goats standing in a pin
Transplantation of sperm-producing stem cells
CRISPR-Cas9 gene editing can improve the effectiveness of spermatogonial stem cell transplantation in mice and livestock, a study finds.
Image credit: Jon M. Oatley.

Similar Articles

Site Logo
Powered by HighWire
  • Submit Manuscript
  • Twitter
  • Facebook
  • RSS Feeds
  • Email Alerts

Articles

  • Current Issue
  • Special Feature Articles – Most Recent
  • List of Issues

PNAS Portals

  • Anthropology
  • Chemistry
  • Classics
  • Front Matter
  • Physics
  • Sustainability Science
  • Teaching Resources

Information

  • Authors
  • Editorial Board
  • Reviewers
  • Librarians
  • Press
  • Site Map
  • PNAS Updates

Feedback    Privacy/Legal

Copyright © 2021 National Academy of Sciences. Online ISSN 1091-6490