Abstract

We present a dynamic causal model that can explain context-dependent changes in neural responses, in the rat barrel cortex, to an electrical whisker stimulation at different frequencies. Neural responses were measured in terms of local field potentials. These were converted into current source density (CSD) data, and the time series of the CSD sink was extracted to provide a time series response train. The model structure consists of three layers (approximating the responses from the brain stem to the thalamus and then the barrel cortex), and the latter two layers contain nonlinearly coupled modules of linear second-order dynamic systems. The interaction of these modules forms a nonlinear regulatory system that determines the temporal structure of the neural response amplitude for the thalamic and cortical layers.

The model is based on the measured population dynamics of neurons rather than the dynamics of a single neuron and was evaluated against CSD data from experiments with varying stimulation frequency (1–40 Hz), random pulse trains, and awake and anesthetized animals. The model parameters obtained by optimization for different physiological conditions (anesthetized or awake) were significantly different.

Following Friston, Mechelli, Turner, and Price (2000), this work is part of a formal mathematical system currently being developed (Zheng et al., 2005) that links stimulation to the blood oxygen level dependent (BOLD) functional magnetic resonance imaging (fMRI) signal through neural activity and hemodynamic variables. The importance of the model described here is that it can be used to invert the hemodynamic measurements of changes in blood flow to estimate the underlying neural activity.

1.  Introduction

The technique of BOLD fMRI is based on the fact that localized neural activity induces co-localized changes in cerebral blood flow (CBF), cerebral metabolic rate of oxygen consumption (CMRO2), and cerebral blood volume (CBV), the set of physiological responses collectively referred to as the hemodynamic response to activation (Kwong et al., 1992; Ogawa et al., 1992). However, while providing important noninvasive spatial information, it is difficult to interpret the blood oxygen level dependent (BOLD) signal in terms of the underlying neural activity. To address this, biophysical models of the chain of processes linking the stimulus and the fMRI signal are being developed.

Among the first to suggest a mathematical forward model of the links between neural activity and the BOLD signal was Friston and colleagues (Friston, Mechelli, Turner, & Price, 2000). Later research (Boas, Jones, Devor, Huppert, & Dale, 2008; Buxton, Uludag, Dubowitz, & Liu, 2004; Huppert, Allen, Benav, Jones, & Boas, 2007; Kong et al., 2004; Martindale et al., 2005; Zheng et al., 2002), added improved revisions and extensions of the component modules that included dynamics. The input to the original Friston system is the neural activity signal, specified directly as a square pulse, thus enabling studies about how aspects of the neural activity affect the observed BOLD signal. However, what this input represents is unclear, and this system leaves the finer details of the relationship between stimulus and the dynamics of neural activity uncharacterized.

A simple nonlinear dynamic model linking stimulus to the amplitude of neural responses was developed by Buxton et al. (2004). The profile of the neural response train was modeled as the difference between an excitatory signal, assumed equivalent to the stimulus, and an inhibitory signal, modeled as a first-order dynamic system with a threshold. However, again, this model is not based on experimental data but mimics the profile of the neural firing rate.

Of particular relevance here is the recent emphasis on the interpretation of event-related experiments that assume both linearity in the hemodynamic response and constancy of the neural response. For this reason, a model capable of representing the temporal dynamics and changes of amplitudes of the neural pulse train is developed and presented in this letter. The model simulates neural responses in the rat barrel cortex to electrical whisker pad stimulation. Its structure is based on the sensory physiological pathway from the brain stem to the thalamus, and then to the cortex, and it may be thought of as a hybrid phenomenological-physiological model of the underlying neural circuitry. The advantage of developing a descriptive model is that the mathematical structure remains simple.

The output of the model is validated against the neural activity data assessed by current source density (CSD) analysis of local field potentials (LFPs) in response to a stimulation pulse train over a range of frequencies (1–40 Hz), the response to random irregular pulse trains, and also two different anesthetic conditions (asesthetized or awake). A single set of parameters can be found to fit well over the above range of stimuli within each anesthetic condition. However, unsurprisingly, across the two different anesthetic states, although the model structure remains invariant, the parameters are significantly different. This context-dependent change in the model parameters may be related to other neural modeling, such as the general dynamic causal model (DCM) framework established by Friston, Harrison, and Penny (2003), where the connections between different areas of the brain are modeled by bilinear parameters reflecting a change in the context.

As the model is expressed in dynamic equations and has relatively few free parameters, it can be easily incorporated into the overall system modeling the link between stimulation, neural activity, and the hemodynamic responses underlying the BOLD signal, enabling us to make inferences about changes in neural activity from the measured BOLD signal. In addition, our model and potential variants also serve as generative models for electrophysiological responses. In this letter, we focus on the LFPs, but it would be possible to augment our model with an electromagnetic forward model to describe electroencephalographic (EEG) signals observed noninvasively at the scalp. There is currently a lot of interest in the use of neural mass models to explain EEG and magnetoencephalographic (MEG) data (David, Harrison, & Friston, 2005; David & Friston, 2003; Garrido, Kilner, Kiebel, Stephan, & Friston, 2007; Jansen & Rit, 1995; Moran et al., 2007), as well as in the neuroimaging context (Almeida & Stetter, 2002; Aubert & Costalat, 2002; Horwitz, Friston, & Taylor, 2000). Unlike established DCMs of EEG and MEG signals, we explicitly model thalamocortical interactions (cf. Robinson, Rennie, Rowe, O'Connor, & Gordon, 2005). However, the model neither describes the exact neural circuitry nor claims to generate neural responses the way the physiological system does. It may nonetheless provide useful insights toward a better understanding of sensory mechanisms such as adaptation as described by Castro-Alamancos (2004b).

2.  Characteristics of the CSD Data

Electrophysiological recordings from Martin, Martindale, Berwick, and Mayhew (2006) have been used for validation of the model. (Readers are directed to their work for details on physiological preparation, stimulus presentation, and data acquisition). The subsequent CSD analysis has been performed as described in Martindale et al. (2003). Although the CSD data from the sink are negative, representing the summed excitatory postsynaptic potential, they were inverted, without loss of generality, so that the amplitude of the CSD response is positive with respect to baseline conditions. Also, the CSD time series were normalized with respect to the magnitude of the first pulse to allow comparison across data sets.

Electrical stimuli were delivered to the whisker pad. For each experimental session under anesthetic condition, 2 second stimuli of 1, 2, 3, 4, 5, 10, 20, and 40 Hz were randomly interleaved. Data were averaged over five animals. For awake experiments, only 1, 2, 5, 10, 20, and 40 Hz stimulus frequencies were used, and data were averaged over six animals. Figure 1 shows the average CSD data for anesthetized animals to trains of different stimulation frequencies. As can be seen, the CSD time series are trains of pulses following the stimulation pulse train, although the amplitude of the individual response varies systematically. At sufficiently high stimulation frequencies, the neural responses to a stimulation pulse train could be characterized in three stages. The neural response to the first pulse of the stimulation train had the highest amplitude, the responses elicited after the first response but before approximately 0.4 s of stimulation were very much reduced in amplitude, and the subsequent responses recovered to a steady state of intermediate amplitude. The amplitude of this steady state was frequency dependent and decreased when the stimulus frequency increased.

Figure 1:

Normalized CSD data in response to stimulations at different frequencies. The steady-state amplitude of the neural responses decreases when the stimulation frequency increases.

Figure 1:

Normalized CSD data in response to stimulations at different frequencies. The steady-state amplitude of the neural responses decreases when the stimulation frequency increases.

Thus, the characteristics of the CSD data across the different stimulation frequencies suggest an underlying adaptation mechanism whose dynamics are much slower than those generating the CSD pulses. The mechanism is activated after each stimulus pulse and subsides within 1 to 2 seconds. Hence for sufficiently low frequencies of stimulation, neural responses are not attenuated. At stimulation frequencies higher than 1 Hz, there is not sufficient time between consecutive stimulus onsets for the mechanism to return to its baseline value. The net effect is a rapid increase in attenuation of neural responses before it decreases to a steady state at a level dependent on the frequency of stimulation.

In addition to regular stimulation pulses, a random stimulus pulse sequence was used to test the model. The stimulus input was constructed by generating a random sequence of integers between 1 and 10, corresponding to frequencies, and the time interval between consecutive pulses is the inverse frequency (e.g., if the integer is 8, then the time interval between the current pulse and the next one is 1/8 second). The duration of the stimulation was 16 s. The averaged CSD data sequence was obtained from 12 anesthetized animals, each of which underwent 30 trials.

3.  Characteristics of the Physiological System

Physiological studies on the whisker-to-barrel pathway (Henderson & Jacquin, 1995; Land, Buffer, & Yaskosky, 1995; Woolsey & Van der Loos, 1970) show that when a stimulus is presented to a rat's whisker, the first area of the brain to respond is the brain stem nuclei (BSN), also called simply brain stem. The brain stem then triggers the thalamus, which sends the information to the cortex.

3.1.  Brain Stem.

Research on the brain stem (Ahissar, Sosnik, & Haidarliu, 2000; Shipley, 1974; Simons, 1985; Sosnik, Haidarliu, & Ahissar, 2001) suggests that the brain stem nuclei act as a relay structure and faithfully transmit the stimulation from the whiskers to the thalamus, with a time delay. In the case of a repetitive stimulation, all responses elicited in the brain stem have the same amplitude, and this amplitude is constant regardless of the stimulation frequency.

3.2.  Thalamus.

However in the thalamus, studies (Castro-Alamancos, 2002b; Chung, Li, & Nelson, 2002; Hartings, Temereanca, & Simons, 2003) have shown that neural responses to a stimulation pulse train decrease over time and converge to a steady state. The amplitude of this steady state decreases as the frequency of stimulation increases. This attenuation of the neural response amplitude in the thalamus can be fitted by an exponential decay (see Figure 2a). The individual shape of neural responses in terms of the histogram of spiking activity can be found in Hartings's work (Hartings et al., 2003), shown as Figure 3a. It can be seen that there is a slight delay in the thalamic response to stimulus onset and that the response lasts for about 15 ms.

Figure 2:

Typical profile of the neural response train in the (a) thalamus and (b) cortex, showing a monotonically decreasing profile in the thalamus but an initial transient characteristic in the cortex.

Figure 2:

Typical profile of the neural response train in the (a) thalamus and (b) cortex, showing a monotonically decreasing profile in the thalamus but an initial transient characteristic in the cortex.

Figure 3:

Typical shape of the neural response to a single stimulation pulse. (a) Figure from Hartings et al. (2003) depicting the histograms of spiking activity in the thalamus in response to whisker deflection at 1, 4, 12, 20, and 40 Hz. Responses to each pulse in a train were averaged, excluding the first. It can be seen that the response lasts around 15 ms. (b) Individual CSD response in the cortex obtained from LFP data in our laboratory from anesthetized animals (black line) and awake animals (gray line). The CSD pulse width is about 20 ms for the anesthetized animals but 10 ms for the awake animals. The subsequent undershoot in the awake animals is also faster and more pronounced than in the anesthetized animals.

Figure 3:

Typical shape of the neural response to a single stimulation pulse. (a) Figure from Hartings et al. (2003) depicting the histograms of spiking activity in the thalamus in response to whisker deflection at 1, 4, 12, 20, and 40 Hz. Responses to each pulse in a train were averaged, excluding the first. It can be seen that the response lasts around 15 ms. (b) Individual CSD response in the cortex obtained from LFP data in our laboratory from anesthetized animals (black line) and awake animals (gray line). The CSD pulse width is about 20 ms for the anesthetized animals but 10 ms for the awake animals. The subsequent undershoot in the awake animals is also faster and more pronounced than in the anesthetized animals.

3.3.  Cortex.

In the cortex, attenuation and undershoot of neural responses in the form of the CSD data are widely observed by researchers (Castro-Alamancos, 2004a; Martindale et al., 2003; Nielsen & Lauritzen, 2001; Sheth et al., 2003). The cortical responses tend to decrease very quickly and then recover to a steady-state amplitude, as shown in Figure 2b. A single CSD cortical response in anesthetized animals is shown in Figure 3b (black line). The response lasts around 20 ms, with a delay of around 5 ms after each stimulus onset (response latency). All responses are followed by an undershoot lasting around 0.2 s. These characteristics seem largely invariant across sessions, animals, or stimulus frequencies (not shown). The cortical neural response of awake animals to an electrical pulse stimulation tends to have faster dynamic characteristics, shown in Figure 3b as a gray line. The response lasts about 10 ms, followed by a more pronounced undershoot.

In what follows, we present a nonlinear dynamic causal model linking the stimulus pulse train to changes in the cortical neural activity represented by the CSD time series. The model will be capable of predicting the characteristics of the neural responses described above.

4.  The Mathematical Model

Like the neural pathway, the model presented here consists of three layers, each corresponding to the three anatomic stations brain stem, thalamus, and cortex respectively, as shown in Figure 4. The stimulus is received by the brain stem, which relays it to the thalamus, which in turn stimulates the cortex.

Figure 4:

Block diagram of the model structure showing three layers corresponding to the brain stem, the thalamus, and the cortex, respectively. The individual thalamic (Th) and cortical (Co) modules are linear second-order systems.

Figure 4:

Block diagram of the model structure showing three layers corresponding to the brain stem, the thalamus, and the cortex, respectively. The individual thalamic (Th) and cortical (Co) modules are linear second-order systems.

The brain stem is modeled as a pure transport delay. Its output represents the average pulse density of action potentials fired by neurons. Similarly, the output from the thalamus to the cortex also represents the average pulse density of action potentials. The difference in amplitude of the output pulses reflects the adaptation characteristics in the thalamus. In contrast to brain stem and thalamus, the output from the cortex models the subthreshold synaptic activity represented by the CSD data in the sink located in layer IV of the cortex.

The mathematical equation used to model the brain stem is the pure transport delay equation:
formula
4.1
where I(t) and zB(t) are the input stimulus and the output from the brain stem, respectively. The time delay is given by τ, that is, the output from the brain stem is the same as its input apart from a pure time delay.

The thalamus is modeled by two linear dynamic modules (Th1 and Th2) and a static nonlinearity fNL. The module Th1 models the shape of the individual pulse density of the spiking activity to a single stimulation pulse. Its output is then modulated in a nonlinear fashion by the output from the second thalamic module Th2, which models the process of sensory adaptation in the feedback loop. It was found that this simple model structure is sufficient to provide a sequence of pulses whose profile decays exponentially.

To capture the characteristics of the CSD data, the model structure of the cortex took on a more complicated form than that of the thalamus, with two additional dynamic modules. Like the thalamus, the module Co1 models the shape of the individual cortical CSD responses, and the feedback module Co2 models the effect of sensory adaptation in the cortex. The additional module Co3 acts as a feedforward inhibition module (Castro-Alamancos, 2004b), which is necessary in order to model the initial transient characteristics observed in the profile of the CSD time series. Finally the undershoot characteristics observed in the cortical CSD data are modeled by the feedforward module Co4. Thus, each module contributes a different dynamic characteristic of the observed neural data rather than a particular physiological process. These modules are coupled (see Figure 4) to form a model relating thalamic neural responses and cortical neural responses.

Each dynamic module is modeled by a linear second-order system, a system that proved to offer a good balance between the number of parameters and the complexity and adaptability of the system response. Linear second-order systems are the main constituent of neural mass models that are widely used in the literature to emulate population dynamics. A nice example of this model is that used by Jansen and Rit (1995) to model oscillations as measured with the EEG. Furthermore, these equations are the same as those used for each population in DCMs of distributed responses. Here we will not focus on the interpretation of the parameters in terms of synaptic rate constants and the like but frame them in terms of their overall impact on the dynamics of a neural population.

The two second-order systems in the thalamus are given by:
formula
4.2
formula
4.3
where eT and sT are the outputs from Th1 and Th2 modules, respectively, with Ki, ζi, and ωi representing the gain, the damping ratio, and the undamped natural frequency for each second-order system. The interaction between the Th1 and Th2 modules is modeled as a nonlinear function, such that the output from Th2 module (sT) attenuates the output from the Th1 module eT as follows:
formula
4.4
where zT is the output signal from the thalamus and is the gain of the Th1 module if sT is zero. The term ( can be regarded as a variable gain of the Th1 module dependent on the amount of sensory adaptation present in the thalamus.
The mathematical equations governing the dynamic cortical modules are similar to those used in the thalamus, with the subscript T replaced by C:
formula
4.5
formula
4.6
formula
4.7
formula
4.8
formula
4.9
formula
4.10
where eC, sC, nC, and uC are the outputs from Co1, Co2, Co3, and Co4 modules, respectively; mC the modulated signal; and zC the model output from the cortical layer.

Each linear second-order system has three unknown parameters (Ki, ζi, and ωi). The overall model has six such systems, yielding 18 parameters to estimate, in addition to the duration of the transport delay. However, with appropriate assumptions and certain invariant properties of the observed neural responses to a single pulse, the modules Th1, Co1, Co3, and Co4 were optimized manually and their parameters fixed once found. The remaining two modules Th2 and Co2 (with six unknown parameters) were then optimized computationally. The optimization procedure is detailed below.

5.  Results

The model was first evaluated against the CSD data collected using regular stimulus pulse trains. It was then tested against CSD data obtained in response to randomly interleaved stimulation pulses.

5.1.  Parameter Adjustments and Model Predictions for Regular Pulse Trains.

The brain stem was modeled by a transport delay of 3 ms under anesthetized conditions and 2 ms under awake conditions. These values were quite robust across frequencies and sessions and were therefore used without further optimization.

5.1.1.  Parameter Tuning of Modules Th1 and Co1.

Both Th1 and Co1 modules provide the shape of the neural response to an individual stimulus pulse, in the thalamus and cortex, respectively. As the shape of the neural pulse from the thalamus is similar to that of the CSD pulse from the cortex, the same model parameters were used for both modules. Furthermore, under a single anesthetic condition (anesthetized or awake), the shape of the CSD pulses seems largely invariant across animals and stimulation frequencies. Across anesthetic conditions, however, the dynamic characteristics of the neural pulse are different, as seen in Figure 3b. Hence, under each anesthetic condition, the six model parameters and for Th1 and Co1 were manually tuned so that the shape of the model predicted CSD pulses matched closely to that of the CSD data. These model parameters were then fixed once satisfactory values were found.

A further constraint on both of these modules is that their impulse response should not produce an undershoot, that is, the model should be overdamped with ζ>1. This is because the undershoot observed in the CSD time series has completely different dynamic characteristics and is modeled separately. The optimized parameters are shown in Table 1 in the Co1 = Th1 row.

Table 1:
Estimated Parameters of the Model.
AnesthetizedAwake
KζωCut-Off Frequency (rad/s)KζωCut-Off Frequency (rad/s)
Th1 1.6 1.1 450 253 1.075 1.1 1800 1011 
Th2a 37.3 2.7 13.0 2.5 30 1.5 15 5.6 
Co1=Th1 1.6 1.1 450 253 1.075 1.1 1800 1011 
Co2a 46.5 1.4 12.7 5.2 80 51.3 
Co3 0.8 0.7 13 13 0.8 0.5 25 38 
Co4 1.8 1.1 12 0.7 30 30 
AnesthetizedAwake
KζωCut-Off Frequency (rad/s)KζωCut-Off Frequency (rad/s)
Th1 1.6 1.1 450 253 1.075 1.1 1800 1011 
Th2a 37.3 2.7 13.0 2.5 30 1.5 15 5.6 
Co1=Th1 1.6 1.1 450 253 1.075 1.1 1800 1011 
Co2a 46.5 1.4 12.7 5.2 80 51.3 
Co3 0.8 0.7 13 13 0.8 0.5 25 38 
Co4 1.8 1.1 12 0.7 30 30 

a The parameters associated with these modules were optimized using the expectation-maximization algorithm.

5.1.2.  Parameter Tuning of Modules Co3 and Co4.

The undershoot characteristic is present in the cortical CSD data but not in the pulse density of the thalamic spiking activity data (see Figure 3); hence, the undershoot is modeled only in the cortex. Two feedforward modules, Co3 and Co4, were used in the cortical layer of the model for this purpose. As the undershoot characteristic under a single anesthetic condition (anesthetized or awake) also seemed to be largely invariant across animals and stimulation frequencies, the six model parameters and for these two modules were also manually tuned and fixed under each anesthetic condition, their values shown in Table 1, in rows Co3 and Co4, respectively. Figure 5 shows an example of the predicted single CSD pulse superimposed with the cortical CSD data for anesthetized and awake animals, respectively. It can be seen that the model predicts the data satisfactorily.

Figure 5:

Model prediction (dashed line) of individual cortical responses superimposed with the measurement (solid line) for (a) anesthetized and (b) awake animals. In both cases, the model fits the peak and undershoot reasonably well.

Figure 5:

Model prediction (dashed line) of individual cortical responses superimposed with the measurement (solid line) for (a) anesthetized and (b) awake animals. In both cases, the model fits the peak and undershoot reasonably well.

5.1.3.  Optimization of Feedback Modules Th2 and Co2.

The feedback modules Th2 and Co2 are largely responsible for modeling the shape of the profile of the response train from the thalamus and the cortex, respectively. The response profiles in the thalamus seem to decrease monotonically to a steady state (Castro-Alamancos, 2002b; Chung et al., 2002; Hartings et al., 2003), whereas in the cortex, the profile drops following the first response before it recovers to a steady state. Thus, the dynamics of the module Co2 in the cortex is likely to be faster than that Th2 in the thalamus.

The profile of the neural responses reflects the sensory adaptation of the responses. It has been shown that characteristics of the sensory adaptation in the thalamus and the cortex vary greatly across stimulation frequencies, as well as across different anesthetic conditions (Castro-Alamancos, 2002b; Chung et al., 2002; Hartings et al., 2003; Martin et al., 2006). Hence, the two feedback modules Th2 and Co2 were optimized, under each anesthetic condition, using the expectation-maximization (EM) algorithm (Friston et al., 2002) in order to obtain the optimal predicted CSD profiles over all stimulus frequencies. This means that the remaining six model parameters ( in the thalamus and in the cortex) associated with the two modules Th2 and Co2 need to be optimized.

The priors for the six model parameters are shown in Table 2. They were found after the initial trial-and-error optimization to ensure numerical convergence of the optimization procedure. We tested the stability of the resulting optimized parameter values by repeating the optimization using different priors and found the model parameters converge to the same optimal set over a wide range of priors.

Table 2:
Priors Associated with the Parameters Optimized Using the Expectation Maximization Algorithm.
AnesthetizedAwake
KζωKζω
Th2 20 20 20 20 
Co2 20 20 20 50 
AnesthetizedAwake
KζωKζω
Th2 20 20 20 20 
Co2 20 20 20 50 

Over the anesthetized conditions, the EM algorithm yielded satisfactory predictions of the CSD data across all stimulation frequencies 1 to about 40 Hz. However, for awake animals, the optimization procedure could not converge satisfactorily, and the resulting predicted responses to high-frequency stimulation (≥10 Hz) were significantly lower than the data at steady state (see Figure 6). Castro-Alamancos (2004b) reported that awake animals showed much less adaptation to repetitive stimuli than anesthetized animals because of enhanced cortico-thalamic feedback but only at high frequencies.

Figure 6:

Profiles of model predictions for awake animals, without using a threshold. For 20 and 40 Hz stimulations, the predicted neural activity profiles (dashed black lines) are significantly lower than the data profiles (solid gray lines) and are outside the error band (representing ±1 standard deviation of the mean of the data).

Figure 6:

Profiles of model predictions for awake animals, without using a threshold. For 20 and 40 Hz stimulations, the predicted neural activity profiles (dashed black lines) are significantly lower than the data profiles (solid gray lines) and are outside the error band (representing ±1 standard deviation of the mean of the data).

This phenomenon could be accounted for by our model by introducing a simple threshold mechanism at the output of the feedback module Th2 for the awake condition. At low frequencies of stimulation, this threshold was not reached and hence had no effect on the model performance. However, at higher stimulation frequencies, the threshold clamped the output from the Th2 module, resulting in less attenuation at the thalamic output, and hence an increase in the CSD amplitude in the cortex. The threshold value of 0.6 was found to be adequate for predicting the CSD responses at high frequencies.

With the threshold mechanism in place, the optimization procedure converged for the awake conditions across the stimulation frequency range 1 to about 40 Hz. The optimized parameters for modules Th2 and Co2 are shown in Table 1, in rows Th2 and Co2, respectively. Also shown in Table 1 are the cut-off frequencies for each module. Because these modules were modeled as second-order systems, they can be thought of as low-pass filters. The cut-off frequency is therefore the bandwidth of such a system and is defined as the frequency at which the magnitude ratio of output to input is 0.7 of (or 3 dB below) its maximum. The MATLAB function bandwidth computes this 3 dB point and therefore was used in our study to compute the cut-off frequencies.

Figure 7 shows the model predicted profiles for anesthetized and awake conditions across all simulation frequencies. The model predictions fit the measured CSD profiles very well, with the same model parameters used within each anesthetic condition for all frequencies.

Figure 7:

Profiles of model predictions to stimulations at different frequencies for anesthetized (left column) and awake (right column) animals. Solid gray lines represent the average profiles of CSD responses over sessions, broken black lines represent the corresponding predictions from the model, and light gray bands represent ±1 standard deviation of the mean of the data.

Figure 7:

Profiles of model predictions to stimulations at different frequencies for anesthetized (left column) and awake (right column) animals. Solid gray lines represent the average profiles of CSD responses over sessions, broken black lines represent the corresponding predictions from the model, and light gray bands represent ±1 standard deviation of the mean of the data.

5.2.  Predictions of Neural Responses to Random Stimulation.

In section 5.1, we described how the model form or structure was optimized using observed responses to regular trains of stimuli. In this section, we establish the validity of our model by showing that the same form can also predict neuronal responses to irregular or random stimuli. In other words, we use the same model form and reoptimized a small number of parameters to show that this form generalizes in its predictive ability.

The experiments (conducted by a different experimenter) used 12 anesthetized animals, and the random stimulus pulse sequence used is shown in Figure 8a, with the averaged CSD response shown in Figure 8b. As the experiments were performed on anesthetized animals, the model parameters were originally kept the same as for regular stimulations of anesthetized animals for all modules. Although the predicted CSD compared reasonably well with the data (not shown), the model prediction was improved by allowing the two feedback modules, Th2 and Co2, to be optimized. These two modules model the sensory adaptation characteristics of the neural responses and hence are crucial in generating profiles of the model output. The six model parameters (optimized using the EM algorithm) associated with these two modules are shown with the corresponding optimized parameters for the regular stimulus condition in Table 3. The model predicted CSD profile is shown in Figure 8c, superimposed with the measured data. The predictions lie well within the error band of ±1 standard deviation.

Figure 8:

(a) Stimulus sequence of randomly interleaved pulses. Intervals between pulses last between 0.1 and 1 s. (b) Average CSD data in response to the random stimulus sequence. (c) Profile of the model predictions (dashed lines) to the random stimulus sequence, compared to the CSD data profile (solid lines). The light gray band corresponds to ±1 standard deviation of the mean of the data.

Figure 8:

(a) Stimulus sequence of randomly interleaved pulses. Intervals between pulses last between 0.1 and 1 s. (b) Average CSD data in response to the random stimulus sequence. (c) Profile of the model predictions (dashed lines) to the random stimulus sequence, compared to the CSD data profile (solid lines). The light gray band corresponds to ±1 standard deviation of the mean of the data.

Table 3:
Estimated Parameters for Regular and Random Stimulation Pulse Trains.
RegularRandom
KζωCut-Off Frequency (rad/s)KζωCut-Off Frequency (rad/s)
Th2 37.3 2.7 13.0 2.5 44.4 3.6 8.9 1.3 
Co2 46.5 1.4 12.7 5.2 47.6 1.7 12.6 4.0 
RegularRandom
KζωCut-Off Frequency (rad/s)KζωCut-Off Frequency (rad/s)
Th2 37.3 2.7 13.0 2.5 44.4 3.6 8.9 1.3 
Co2 46.5 1.4 12.7 5.2 47.6 1.7 12.6 4.0 

5.3.  Results Summary.

The model has been shown to be capable of capturing the temporal characteristics of the individual and amplitude profiles of CSD data for both anesthetized and awake animals and for regular and random stimuli. As 12 of the 18 model parameters can be fixed, fitting new data requires only optimizing the remaining six parameters associated with the feedback modules Th2 and Co2.

6.  Discussion

A nonlinear dynamic causal model that relates changes in neural activity (in terms of the CSD of LFP) to stimulation pulse trains over a range of frequencies has been evaluated. This section examines the structure of the model and provides some physiological interpretations of its modules where possible. The generalization of the model to other stimulus types and its limitations are also discussed.

6.1.  Structure of the Model.

The first organizational layer of the model is structured as the whisker-to-barrel neural pathway with three stations: brain stem, thalamus, and cortex. The brain stem is treated as a simple relay with delay. Thus, the input to the thalamus is simply a delayed version of the input pulse train.

The thalamus is modeled by the nonlinear interaction of two dynamic modules, Th1 and Th2. The output of Th1, modulated by the output of Th2, may be considered to be the average pulse density of action potentials transmitted to the cortex. Again this is represented as an essentially smoothed pulse train but with modulated amplitude. The module Th1 may be regarded as modeling the excitation component of the neural responses, with module Th2 providing neural feedback, suppression, or modulation as observed in the data. In the absence of Th2, the thalamic response to a stimulus pulse train of uniform amplitude would show no amplitude modulation. It is worthwhile noting that the neural suppression may be due to various causes. In the thalamus, the depression of synapses, the inhibitory action of the reticular nuclei, and the feedback from the cortex may be involved (Castro-Alamancos, 2004b; Chung et al., 2002; Hartings et al., 2003). These sources, however, are not modeled separately by the current model. The overall effect of the “suppression” is modeled by the single module Th2. The dynamic characteristics of Th2 are much slower than those of Th1. This is reflected in the very different cut-off frequencies of the two modules shown in Table 1.

In the cortex, further suppression (or attenuation of the amplitude profile) is observed and may be due to internal inhibitory networks (Beierlein, Gibson, & Connors, 2003). Some of these mechanisms are part of a phenomenon termed sensory adaptation, which is mainly present during anesthetized states and results in a reduction over time in the responsiveness of a sensory system to a periodic stimulation (Castro-Alamancos, 1997, 2004a). The overall effect of this suppression on the amplitude profile, whatever the mechanisms involved, is modeled by the feedback module Co2.

A comparison of the model predicted profiles of thalamic and cortical neural responses, shown in Figure 9, suggests that most of the attenuation observed in the cortical CSD signal at steady state already exists at the level of the thalamic output.

Figure 9:

Model thalamic (dashed line) and cortical (solid line) outputs for a 5 Hz stimulation in the anesthetized case, demonstrating that most of the attenuation observed at steady state in the cortical output is already present in the thalamic output.

Figure 9:

Model thalamic (dashed line) and cortical (solid line) outputs for a 5 Hz stimulation in the anesthetized case, demonstrating that most of the attenuation observed at steady state in the cortical output is already present in the thalamic output.

An additional feature present in the cortical CSD data is the poststimulus undershoot in the shape of the individual neural response. The undershoot may represent the hyperpolarization of neuron cell bodies following the arrival of synchronized excitatory postsynaptic potentials (EPSPs). Although this feature is modeled by two feedforward modules Co3 and Co4, module Co3 was found to also play a crucial part in modeling the profile shape of the CSD data because its output partly drives the feedback module Co2, which models the suppression in the cortex. Finally, as in the model of the thalamus, the cortical module Co1 can be regarded as modeling the excitatory component of the neural responses.

6.2.  Awake Versus Anesthetized States.

An important difference between anesthetized and awake data is the degree of attenuation at steady state. In the awake condition, the steady-state amplitude is much higher across all stimulation frequencies. An explanation of this was given by Castro-Alamancos (2004b): awake animals have a high baseline of adaptation, consequently are less capable of further adaptation to repetitive stimulation. In contrast, a greater degree of suppression at steady state is observed for anesthetized animals, for which the adaptation baseline level is low.

The spike and the undershoot characteristics of the CSD data are also markedly different for the two anesthetic conditions. In the awake condition, the spike and the undershoot are faster than for anesthetized animals, implying different mechanisms at the individual cell level and increased synchrony in neuron firing for awake animals.

The above characteristics of the data are reflected in the optimized model parameters, resulting in higher cut-off frequencies for all modules under the awake condition, shown in Table 1. This can be interpreted as a result of a dampening effect on neural activity during anesthesia. As shown in many studies (Berwick et al., 2002; Lahti, Ferris, Li, Sotak, & King, 1999; Martin et al., 2002, 2006; Peeters, Tindemans, De Schutter, & Van der Linden, 2001; Shtoyerman, Arieli, Slovin, Vanzetta, & Grinvald, 2000; Sicard et al., 2003), this dampening effect of anesthesia is also observed in the CBF and the BOLD signals.

6.3.  Suppression Threshold in Awake Animals.

According to Castro-Alamancos (2004b) in the awake animal, additional mechanisms come into effect when the frequency of stimulation is sufficiently high. As described above, we found it necessary to clamp the output of the thalamic feedback module in order to fit the data at high frequencies (>10 Hz). The model predicted steady-state amplitude of the response train to high-frequency stimulations is much improved by this modification. Castro-Alamancos (2002a, 2002b) and Fanselow and Nicolelis (1999) showed that when the animal is awake, high-frequency signals are better transmitted, suggesting the activation of an unspecified facilitation mechanism. Facilitating the transmission of the stimulus when it reaches high frequencies is an emergent effect of clamping modulatory inhibition and may reflect the underlying physiological process.

6.4.  Impact of the Stimulation Pulse Width.

All electrophysiological experiments presented in this letter were conducted with a stimulus pulse width of 0.3 ms. However, to reduce the computational load, the input pulse train used in our model has a pulse width of 5 ms. Studies have shown that the neural response width is invariant to the stimulus pulse width, provided it is fewer than a few milliseconds. Rosengarten, Lutz, and Hossmann (2003), using stimulus pulse widths of 0.3 ms and 5 ms, observed that the width of the evoked field potentials was not affected. Similar studies were carried out in our laboratory (Devonshire, unpublished results) in which stimulation pulse widths of 0.3, 1, 3, 5, and 6.5 ms were tested, and CSD and optical imaging spectroscopy (OIS) data were collected concurrently. The results demonstrated that the neural responses do not change with the stimulus pulse width within the range tested, in agreement with the observation by Rosengarten et al. (2003).

That study and the CSD data presented in this letter suggest that the use of 5 ms pulses in the input pulse train of the model is a justifiable approximation for all short-duration pulses.

6.5.  Response to Pulse Stimuli of Various Intensities.

For the data used in this letter, the anesthetized and awake animals were presented with stimuli of intensity of 1.2 mA and 0.3 mA, respectively. All data were normalized in order to compare different sessions and animals, and the model was developed to fit these normalized data in response to a unitary input pulse train. For a linear dynamic system, the stimulus intensity does not affect the model parameters, and hence it is normal practice to use normalized variables. However, our model is nonlinear, and the implication is that the model may not be generalized for different stimulus intensities. This was investigated by testing the model on several CSD data sets obtained under stimulus intensities ranging from 0.6 mA to 1.6 mA (Hewson-Stoate, Jones, Martindale, Berwick, & Mayhew, 2005; Jones, Hewson-Stoate, Martindale, Redgrave, & Mayhew, 2004) for anesthetized animals. The model was found to be able to predict the CSD data well at different stimulus intensities but with different values of the model parameters. Thus, although the model structure seems invariant over stimulus intensities, a nonlinear scaling may be required to model responses to different stimulation intensities.

7.  Conclusion and Future Development

The dynamic causal model described in this letter is a nonlinear lumped-parameter model that has a three-station structure based on the physiology of the whisker-to-barrel pathway and incorporates some physiological aspects of their behavior. The construction of the model follows neither a systematic physiological nor a systematic signal analysis approach. We are not modeling individual neurons; neither are we modeling a generic black box 18 parameter filter. The overall structure of the model is guided by the physiological process, but each subsystem is a lumped-parameter model, capturing a particular dynamic characteristic in the neural signal. Some of the parameters are used to model the temporal characteristics of the summed depolarization and the hyperpolarization of a population of neurons, while others are involved in the modulation of the neural response profiles by suppression and adaptation. The modular structure of the model, and its decomposition into the different subsystems, provides important restrictions on the possible interactions between the model parameters during the optimization procedure. The advantage of such a modeling approach is that it reduces the complexity of the model and allows physical interpretation of the model parameters.

Under different anesthetic regimes, components of the model can be reoptimized while the structure of the model remains unchanged. The model can thus be used to simulate neural responses under anesthetized or awake conditions, over a range of stimulation frequencies, and to regular or random stimuli.

Several forward models have been developed to relate neural responses to the BOLD signal through changes in cerebral blood flow (CBF), cerebral blood volume (CBV), oxygen extraction fraction (OEF), and deoxyhemoglobin (Hbr). The ultimate aim of this modeling effort is to estimate the evoked changes in neural activity from the measured BOLD signal. Our model linking stimulus to neural activity can be incorporated into such forward models. A method to use the model effectively in the context of fMRI interpretation will be presented in a future article.

References

Ahissar
,
E.
,
Sosnik
,
R.
, &
Haidarliu
,
S.
(
2000
).
Transformation from temporal to rate coding in a somatosensory thalamocortical pathway
.
Nature
,
406
(
6793
),
302
306
.
Almeida
,
R.
, &
Stetter
,
M.
(
2002
).
Modeling the link between functional imaging and neuronal activity: Synaptic metabolic demand and spike rates
.
Neuroimage
,
17
(
2
),
1065
1079
.
Aubert
,
A.
, &
Costalat
,
R.
(
2002
).
A model of the coupling between brain electrical activity, metabolism, and hemodynamics: Application to the interpretation of functional neuroimaging
.
Neuroimage
,
17
(
3
),
1162
1181
.
Beierlein
,
M.
,
Gibson
,
J. R.
, &
Connors
,
B. W.
(
2003
).
Two dynamically distinct inhibitory networks in layer 4 of the neocortex
.
J. Neurophysiol.
,
90
(
5
),
2987
3000
.
Berwick
,
J.
,
Martin
,
C.
,
Martindale
,
J.
,
Jones
,
M.
,
Johnston
,
D.
,
Zheng
,
Y.
, et al
(
2002
).
Hemodynamic response in the unanesthetized rat: Intrinsic optical imaging and spectroscopy of the barrel cortex
.
J. Cereb. Blood Flow Metab.
,
22
(
6
),
670
679
.
Boas
,
D. A.
,
Jones
,
S. R.
,
Devor
,
A.
,
Huppert
,
T. J.
, &
Dale
,
A. M.
(
2008
).
A vascular anatomical network model of the spatio-temporal response to brain activation
.
Neuroimage
,
40
(
3
),
1116
1129
.
Buxton
,
R. B.
,
Uludag
,
K.
,
Dubowitz
,
D. J.
, &
Liu
,
T. T.
(
2004
).
Modeling the hemodynamic response to brain activation
.
Neuroimage
,
23
(
Suppl. 1
),
S220
233
.
Castro-Alamancos
,
M. A.
(
1997
).
Short-term plasticity in thalamocortical pathways: Cellular mechanisms and functional roles
.
Rev. Neurosci.
,
8
(
2
),
95
116
.
Castro-Alamancos
,
M. A.
(
2002a
).
Different temporal processing of sensory inputs in the rat thalamus during quiescent and information processing states in vivo
.
J. Physiol.
,
539
(
Pt. 2
),
567
578
.
Castro-Alamancos
,
M. A.
(
2002b
).
Properties of primary sensory (lemniscal) synapses in the ventrobasal thalamus and the relay of high-frequency sensory inputs
.
J. Neurophysiol.
,
87
(
2
),
946
953
.
Castro-Alamancos
,
M. A.
(
2004a
).
Absence of rapid sensory adaptation in neocortex during information processing states
.
Neuron
,
41
(
3
),
455
464
.
Castro-Alamancos
,
M. A.
(
2004b
).
Dynamics of sensory thalamocortical synaptic networks during information processing states
.
Prog. Neurobiol.
,
74
(
4
),
213
247
.
Chung
,
S.
,
Li
,
X.
, &
Nelson
,
S. B.
(
2002
).
Short-term depression at thalamocortical synapses contributes to rapid adaptation of cortical sensory responses in vivo
.
Neuron
,
34
(
3
),
437
446
.
David
,
O.
, &
Friston
,
K. J.
(
2003
).
A neural mass model for MEG/EEG: Coupling and neuronal dynamics
.
Neuroimage
,
20
(
3
),
1743
1755
.
David
,
O.
,
Harrison
,
L.
, &
Friston
,
K. J.
(
2005
).
Modeling event-related responses in the brain
.
Neuroimage
,
25
(
3
),
756
770
.
Fanselow
,
E. E.
, &
Nicolelis
,
M. A.
(
1999
).
Behavioral modulation of tactile responses in the rat somatosensory system
.
J. Neurosci.
,
19
(
17
),
7603
7616
.
Friston
,
K. J.
,
Harrison
,
L.
, &
Penny
,
W.
(
2003
).
Dynamic causal modeling
.
Neuroimage
,
19
(
4
),
1273
1302
.
Friston
,
K. J.
,
Mechelli
,
A.
,
Turner
,
R.
, &
Price
,
C. J.
(
2000
).
Nonlinear responses in fMRI: The balloon model, Volterra kernels, and other hemodynamics
.
Neuroimage
,
12
(
4
),
466
477
.
Friston
,
K. J.
,
Penny
,
W.
,
Phillips
,
C.
,
Kiebel
,
S.
,
Hinton
,
G.
, &
Ashburner
,
J.
(
2002
).
Classical and Bayesian inference in neuroimaging: Theory
.
Neuroimage
,
16
(
2
),
465
483
.
Garrido
,
M. L.
,
Kilner
,
J. M.
,
Kiebel
,
S. J.
,
Stephan
,
K. E.
, &
Friston
,
K. J.
(
2007
).
Dynamic causal modeling of evoked potentials: A reproducibility study
.
Neuroimage
,
36
(
3
),
571
580
.
Hartings
,
J. A.
,
Temereanca
,
S.
, &
Simons
,
D. J.
(
2003
).
Processing of periodic whisker deflections by neurons in the ventroposterior medial and thalamic reticular nuclei
.
J. Neurophysiol.
,
90
(
5
),
3087
3094
.
Henderson
,
T. A.
, &
Jacquin
,
M. F.
(
1995
).
What makes subcortical barrels? Requisite trigeminal circuitry and developmental mechanisms
. In
E. G. Jones & I.T. Diamond (Eds.)
,
The barrel cortex of rodents
(pp. 123–187)
.
New York
:
Plenum Press
.
Hewson-Stoate
,
N.
,
Jones
,
M.
,
Martindale
,
J.
,
Berwick
,
J.
, &
Mayhew
,
J.
(
2005
).
Further nonlinearities in neurovascular coupling in rodent barrel cortex
.
Neuroimage
,
24
(
2
),
565
574
.
Horwitz
,
B.
,
Friston
,
K. J.
, &
Taylor
,
J. G.
(
2000
).
Neural modeling and functional brain imaging: An overview
.
Neural Netw.
,
13
(
8–9
),
829
846
.
Huppert
,
T. J.
,
Allen
,
M. S.
,
Benav
,
H.
,
Jones
,
P. B.
, &
Boas
,
D. A.
(
2007
).
A multicompartment vascular model for inferring baseline and functional changes in cerebral oxygen metabolism and arterial dilation
.
J. Cereb. Blood Flow Metab.
,
27
(
6
),
1262
1279
.
Jansen
,
B. H.
, &
Rit
,
V. G.
(
1995
).
Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns
.
Biol. Cybern.
,
73
(
4
),
357
366
.
Jones
,
M.
,
Hewson-Stoate
,
N.
,
Martindale
,
J.
,
Redgrave
,
P.
, &
Mayhew
,
J.
(
2004
).
Nonlinear coupling of neural activity and CBF in rodent barrel cortex
.
Neuroimage
,
22
(
2
),
956
965
.
Kong
,
Y.
,
Zheng
,
Y.
,
Johnston
,
D.
,
Martindale
,
J.
,
Jones
,
M.
,
Billings
,
S.
, et al
(
2004
).
A model of the dynamic relationship between blood flow and volume changes during brain activation
.
J. Cereb. Blood Flow Metab.
,
24
(
12
),
1382
1392
.
Kwong
,
K. K.
,
Belliveau
,
J. W.
,
Chesler
,
D. A.
,
Goldberg
,
I. E.
,
Weisskoff
,
R. M.
,
Poncelet
,
B. P.
, et al
(
1992
).
Dynamic magnetic resonance imaging of human brain activity during primary sensory stimulation
.
Proc. Natl. Acad. Sci. U.S.A.
,
89
(
12
),
5675
5679
.
Lahti
,
K. M.
,
Ferris
,
C. F.
,
Li
,
F.
,
Sotak
,
C. H.
, &
King
,
J. A.
(
1999
).
Comparison of evoked cortical activity in conscious and propofol-anesthetized rats using functional MRI
.
Magn. Reson. Med.
,
41
(
2
),
412
416
.
Land
,
P. W.
,
Buffer
,
S. A.
, Jr
, &
Yaskosky
,
J. D.
(
1995
).
Barreloids in adult rat thalamus: Three-dimensional architecture and relationship to somatosensory cortical barrels
.
J. Comp. Neurol.
,
355
(
4
),
573
588
.
Martin
,
C.
,
Berwick
,
J.
,
Johnston
,
D.
,
Zheng
,
Y.
,
Martindale
,
J.
,
Port
,
M.
, et al
(
2002
).
Optical imaging spectroscopy in the unanaesthetized rat
.
J. Neurosci. Methods.
,
120
(
1
),
25
34
.
Martin
,
C.
,
Martindale
,
J.
,
Berwick
,
J.
, &
Mayhew
,
J.
(
2006
).
Investigating neural-hemodynamic coupling and the hemodynamic response function in the awake rat
.
Neuroimage
,
32
(
1
),
33
48
.
Martindale
,
J.
,
Berwick
,
J.
,
Martin
,
C.
,
Kong
,
Y.
,
Zheng
,
Y.
, &
Mayhew
,
J.
(
2005
).
Long duration stimuli and nonlinearities in the neural-haemodynamic coupling
.
J. Cereb. Blood Flow Metab.
,
25
(
5
),
651
661
.
Martindale
,
J.
,
Mayhew
,
J.
,
Berwick
,
J.
,
Jones
,
M.
,
Martin
,
C.
,
Johnston
,
D.
, et al
(
2003
).
The hemodynamic impulse response to a single neural event
.
J. Cereb. Blood Flow Metab.
,
23
(
5
),
546
555
.
Moran
,
R.
,
Kiebel
,
S. J.
,
Stephan
,
S. J.
,
Reilly
,
R. B.
,
Daunizeau
,
J.
, &
Friston
,
K. J.
(
2007
).
A neural mass model of spectral responses in electrophysiology
.
NeuroImage
,
37
(
3
),
706
720
.
Nielsen
,
A. N.
, &
Lauritzen
,
M.
(
2001
).
Coupling and uncoupling of activity-dependent increases of neuronal activity and blood flow in rat somatosensory cortex
.
Journal of Physiology–London
,
533
(
3
),
773
785
.
Ogawa
,
S.
,
Tank
,
D. W.
,
Menon
,
R.
,
Ellermann
,
J. M.
,
Kim
,
S. G.
,
Merkle
,
H.
, et al
(
1992
).
Intrinsic signal changes accompanying sensory stimulation: Functional brain mapping with magnetic resonance imaging
.
Proc. Natl. Acad. Sci. U.S.A.
,
89
(
13
),
5951
5955
.
Peeters
,
R. R.
,
Tindemans
,
I.
,
De Schutter
,
E.
, &
Van der Linden
,
A.
(
2001
).
Comparing BOLD fMRI signal changes in the awake and anesthetized rat during electrical forepaw stimulation
.
Magn. Reson. Imaging.
,
19
(
6
),
821
826
.
Robinson
,
P. A.
,
Rennie
,
C. J.
,
Rowe
,
D. L.
,
O'Connor
,
S. C.
, &
Gordon
,
E.
(
2005
).
Multiscale brain modeling
.
Phil. Trans. R. Soc. B
,
360
,
1043
1050
.
Rosengarten
,
B.
,
Lutz
,
H.
, &
Hossmann
,
K. A.
(
2003
).
A control system approach for evaluating somatosensory activation by laser: Doppler flowmetry in the rat cortex
.
J. Neurosci. Methods
,
130
(
1
),
75
81
.
Sheth
,
S.
,
Nemoto
,
M.
,
Guiou
,
M.
,
Walker
,
M.
,
Pouratian
,
N.
, &
Toga
,
A. W.
(
2003
).
Evaluation of coupling between optical intrinsic signals and neuronal activity in rat somatosensory cortex
.
Neuroimage
,
19
(
3
),
884
894
.
Shipley
,
M. T.
(
1974
).
Response characteristics of single units in the rat's trigeminal nuclei to vibrissa displacements
.
J. Neurophysiol.
,
37
(
1
),
73
90
.
Shtoyerman
,
E.
,
Arieli
,
A.
,
Slovin
,
H.
,
Vanzetta
,
I.
, &
Grinvald
,
A.
(
2000
).
Long-term optical imaging and spectroscopy reveal mechanisms underlying the intrinsic signal and stability of cortical maps in V1 of behaving monkeys
.
J. Neurosci.
,
20
(
21
),
8111
8121
.
Sicard
,
K.
,
Shen
,
Q.
,
Brevard
,
M. E.
,
Sullivan
,
R.
,
Ferris
,
C. F.
,
King
,
J. A.
, et al
(
2003
).
Regional cerebral blood flow and BOLD responses in conscious and anesthetized rats under basal and hypercapnic conditions: Implications for functional MRI studies
.
J. Cereb. Blood Flow Metab.
,
23
(
4
),
472
481
.
Simons
,
D. J.
(
1985
).
Temporal and spatial integration in the rat SI vibrissa cortex
.
J. Neurophysiol.
,
54
(
3
),
615
635
.
Sosnik
,
R.
,
Haidarliu
,
S.
, &
Ahissar
,
E.
(
2001
).
Temporal frequency of whisker movement. I. Representations in brain stem and thalamus
.
J. Neurophysiol.
,
86
(
1
),
339
353
.
Woolsey
,
T. A.
, &
Van der Loos
,
H.
(
1970
).
The structural organization of layer IV in the somatosensory region (SI) of mouse cerebral cortex: The description of a cortical field composed of discrete cytoarchitectonic units
.
Brain. Res.
,
17
(
2
),
205
242
.
Zheng
,
Y.
,
Johnston
,
D.
,
Berwick
,
J.
,
Chen
,
D.
,
Billings
,
S.
, &
Mayhew
,
J.
(
2005
).
A three-compartment model of the hemodynamic response and oxygen delivery to brain
.
Neuroimage
,
28
(
4
),
925
939
.
Zheng
,
Y.
,
Martindale
,
J.
,
Johnston
,
D.
,
Jones
,
M.
,
Berwick
,
J.
, &
Mayhew
,
J.
(
2002
).
A model of the hemodynamic response and oxygen delivery to brain
.
Neuroimage
,
16
(
3 Pt 1
),
617
637
.