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

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 (A_opt,T_opt) 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.

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.