A few previous studies have made estimates of the effects of different neural models on the form of the BOLD response dynamics. For example, Mechelli et al. (2001) report a simulation study of the estimated regional cerebral blood flow (rCBF) and BOLD signals as a function of the duration, onset asynchrony, and relative amplitudes of two brief stimuli. They included a basic model of neuronal dynamics and varied one parameter of this model – the amplitude of a slow late transient – to show its effect on the simulated BOLD responses. This exercise constitutes an unvalidated forward model of the effect of only one parameter of a simulated neural response. However, what they offer as the analysis of the effects on the BOLD waveform is described (incomprehensibly) as “the BOLD parameter estimates,” which is actually a single-valued function with no specification of which parameter(s) is/are being estimated. Since no comparison is made with empirical data or their noise limitations, the fact that this study includes a proposed model for the neural dynamics does not qualify it as a validated procedure for estimating the neural population dynamics underlying the local BOLD signals (which is the goal of our study).

Buxton et al. (2004) extend their balloon model of the hemodynamic response leading to the BOLD signal by proposing a model of the neural response to account for the temporal non-linearity in BOLD responses as a function of duration. This neural response incorporates a slow subtractive inhibitory component to the net neural signal, which has the effect of producing a neural response consisting of an initial transient followed by a sustained plateau. Model responses for three kinds of stimuli – a single short pulse, two short pulses, and one long pulse, are offered as a demonstration of the properties of this model. As with Mechelli et al. (2001), no attempt is made to compare the model outputs with actual BOLD recordings, so the Buxton et al. (2004) study again does not qualify as a validated procedure for estimating the neural population dynamics underlying the local BOLD signals.

Both Mechelli et al. (2001) and Buxton et al. (2004) include parameters intended to account for the temporal non-linearity of short-duration responses (which do not fall linearly as response duration is reduced; Boynton et al., 1996; Birn et al., 2001). Both studies demonstrate the required lack of reduction in a single example of a short-duration response, but neither study provides a validation of either the waveform or the amplitude response function relative to empirical BOLD data. In principle, either model could provide a platform for such validation or the further estimation of the dynamics of the underlying neural population response, but they neither do so nor suggest procedures by which such estimation could be achieved.

Quantitative estimation of the temporal transients underlying BOLD responses can be achieved by the use of an appropriate set of temporal duration stimuli, such as the temporal summation study reported in Birn et al. (2001). Their stimuli consisted of counterphasing checkerboard patches with durations increasing in factors of 2 from 0.5–2 s (together with a 20-s duration stimulus). The key result was that the BOLD response amplitudes did not increase proportionately with stimulus duration. (A proportionate increase would be expected under a general linear model (GLM) of linear summation with a neural response of constant amplitude.) Instead, the amplitude showed little increase from the short durations, implying either a non-linear summation process or a disproportionate strength of neural signals for short-duration stimuli (Figure 2).

Bandettini and Ungerleider (2001) proposed a discrete approximation analysis for the derivation of the implied neural signal in human cortex under the assumption that the integration process underlying the generation of the BOLD response was indeed equivalent to a linear convolution operation. To derive a crude estimate of the neural signal, they took the difference between integrated BOLD response areas for each successive pair of stimulus durations as an estimate of the amplitude of the neural signal for the time interval defined by the difference in stimulus durations (Figure 2), based on the (unstated) assumption that the neural response is fast enough to follow the stimulus exactly. This discrete approximation approach implied that the neural signal underlying these BOLD responses has a pronounced early transient followed by a sustained phase of stable amplitude out to 20 s (red box functions, Figure 2). The details of this transient, however, remained unresolved.

In this section, we propose a more biophysically based development of the forward modeling approach that we term NDFO. Rather than simply characterizing the behavior of the BOLD waveform (Birn et al., 2001; Bandettini et al., 2002; d'Avossa et al., 2003; Fox et al., 2005; Grotz et al., 2009) or attempting to infer the potentially complex properties of the underlying neural mechanisms from the form of the BOLD response by deconvolution (Glover, 1999; Logothetis, 2002, 2003; Logothetis and Wandell, 2004), the concept of forward modeling is to start from the “cause” – the neural signal, instead from the “consequence” – the BOLD signal. Thus we incorporate into the model as much knowledge as possible about the likely neural substrate and then optimize the remaining details to best fit the BOLD waveform. This knowledge includes the known temporal properties of neural responses and non-linearities both of the population of neural responses to the stimulus and of the metabolic requirements of the neural processing that lead to the measurable BOLD response.

The dynamic forward modeling discussed so far is linear (see Boynton et al., 1996) in the sense that it assumes no non-linearity in the model linking the neural response to the BOLD response (or indeed in linking the original stimulus waveform to the BOLD response). It merely estimates the gain of the neural response for each stimulus duration, which can vary by a linear process. A more general form of forward modeling is to incorporate a variety of possible non-linearities into the structure of the model. The range of possible non-linearities at each stage of the model is so large that the approach is not feasible unless one incorporates Bayesian constraints in the modeling, based on the known biophysics of the neuronal response properties likely to underlie the BOLD activation (Logothetis and Wandell, 2004; see Table 1). The particular form of the model will depend on the specific knowledge of the population dynamics of the signal driving the hemodynamic response available from the literature at any given time.

The neural model that we investigate in the present version of the analysis is the sum of positive and a negative component based on delayed gamma functions convolved with the stimulus waveform, with the parametrization specified in Table 2. To illustrate the properties of the model, we analyze the effect of varying the inhibitory ratio implied by the negative component weight b, and the offset/onset gain ratio λ.

The starting point for the NDFO modeling is a model of the neural signal, whose first effect in terms of the cascade of BOLD dynamics is to create a metabolic demand G(t) in the neighboring glial cells (see Figure 1). Gamma probability density functions have the analytic form a·t^k−1·e^−t/τ/Γ(k), where t is the time dimension, a is a scaling parameter, and k and τ are generic waveform parameters. (To equate the area to unity, the function is defined with a = 1/τ^k), but for the present purposes a is a free optimization parameter.) They may be termed simply “gamma functions” to emphasize their analytic rather than statistical properties. In engineering the same function is known as the n-pole filter function (with n = k), and is used to describe the dynamics of a wide range of processes. For the present application, the temporal gamma function is assumed to have integer powers of k and corresponds to the solution of differential equations with real (non-imaginary) roots. The gamma function has the analytic advantage over many other functions, such as the Gaussian, that it is by definition causal because it has the value of 0 at t = 0, and is defined as 0 for t < 0 (i.e., the full specification is y = k^t−1·e^−t/τ for t < 0; y = 0 for t < 0). Its shape progresses from highly asymmetric around the peak for small k to approximately Gaussian and symmetric for large k (see Figure 3).

A key feature of this formalism is that the peak latency is determined by the time constant, τ, and is proportional to the width (at half height) of the response peak, which can be estimated to good accuracy by the methods described in the next section.

A comprehensive model of the BOLD therefore requires an accurate model of the intracellular potential dynamics deriving from the sensory stimulation. Formally, we propose to use the full neural response model jointly specified in Eqs 1 and 2:

where ε is the net source of additive noise and the function is governed by the parameter set.

The neural signal N_P(t) in Eq. 1 is the sum of positive and a negative component based on delayed gamma functions convolved with the stimulus waveform:

where N(t)=(1−t )⋅tkn−1e−t/τn−bτntkne−t/τn

(Note the convention that time series functions are capitalized, impulse response kernels are lower case, and vectors are bold-face.) The neural impulse response expression n(t) is set up so that its convolution with a step function is equivalent to the sum of a pure transient and a pure sustained component. In addition, the expression is specified with an additional transmission delay, Δt, that delays the response relative to the stimulus without affecting its waveshape. Thus, Table 2 defines the parameters in vector p of this equation.

One issue that arises is how to measure the latency Δt of the delayed gamma functions of Figure 3. A simple derivation can show that the peak latency t_peak of these responses is specified by the expression (kτ_n + Δt). Thus if k and Δt, both of which are well-determined from neurophysiological studies in monkey cortex, are set to means of their Bayesian priors on this basis, the peak latency t_peak is determined from the value of τ, which can be accurately derived from the model optimization.

Since little is known about the glial dynamics of transmitter recovery, we will pursue two options as to their effects. One option is that the metabolic demand driving the BOLD response derives from the transmitter recovery cycle following the activation by an axonal spike. Since axonal spikes represent only the positive aspect of the intracellular voltage and since 90% of cortical synapses are excitatory (Shank and Aprison, 1979; Wang and Floor, 1994), the signal transmitted from one cortical stage to the next may be treated as a half-wave rectified version of the dynamic neural signal. This prediction is shown as the blue curves in Figure 4Aa (indexed in row/column notation), which is an overlay of the model estimates of the neural responses to stimulus pulses that double in duration from 8 ms to 16 s (eight doublings). For this example, the neural response has balanced excitation and inhibition, so even the prolonged pulses generate only an initial transient response, with the negative lobe at offset being thresholded out by the rectification. (Note that the local metabolic demand, M_i(t), has the same time course in this model as the transmitter recovery from which it derives. The energetic processes required for the recovery to the initial state, however, form an oxygen-based chain of glial metabolic response, G(t), that may have substantially slower time course at one or more stages.)

The other option is to consider the instantaneous metabolic demand of both excitatory and inhibitory cells, or both “on” and “off” cells, implying that the signal generating the metabolic demand is a full-wave rectified version of the intracellular voltage. This option is shown as the red curve in Figure 4Aa, whose spikes represent the instantaneous metabolic demand for the offset transient for each of the pulse duration doublings (which were thresholded out in the blue curve).

Having provided relevant variables for the linear and non-linear components of the neural response dynamics, we may now consider the issue of the metabolic coupling of the neural signal to generate the BOLD response. This coupling has been extensively modeled over the past two decades, with the best-known example being the Buxton–Friston balloon model and the most elaborated version being by Sotero and Trujillo-Barreto (2007, 2008). However, these models are not well-validated by empirical human studies (because they contain too many interdependent variables to allow the independent assessment of each separately.

We may now evaluate the response to these two options for the non-linearity of the metabolic demand through the biophysical chain of the metabolic processes to the measured BOLD signal (Table 1). The first element in this chain is the astrocytes surrounding the neuron, which provide glucose to the neuron and replenish its supply by ATP metabolism fueled by oxygen from neighboring blood vessels. It may be emphasized that the astrocyte metabolic processes are as slow, relative to the intracellular signal dynamics, as are the processes of hemodynamic oxygen supply. The time constant of the astrocyte responses at the cell body is known to be of the order of several seconds (Kelly and van Essen, 1974; Filosa et al., 2004; Metea and Newman, 2006), and it is clear that these slow responses must mediate the hemodynamic component. However, at present too little is known of the dynamics of transmitter recovery and/or the non-linearities in the process to securely assign time constants to the astrocytic component relative to the hemodynamic component of the metabolic coupling. For the present example, we will therefore treat the entire chain from the metabolic demand to the magnetic resonance signal in the traditional fashion, as a unitary linear kernel. As stated above, this kernel is often termed the HRF, but in view of its likely substantial astrocyte contribution, we give it the more general term of the MRK. Our main goal is to estimate the properties of the neural signal processing, and it will be seen that there is sufficient information to provide a rich analysis of these properties, and to account for the empirical non-linearities of the BOLD signal, as long as the metabolic supply chain conforms to the linearity assumption. (As more information becomes available, non-linear aspects of the metabolic coupling may readily be incorporated into the analysis we propose.)

The development of the peak amplitude in the BOLD temporal summation series is shown as the blue summation curve in Figure 4Ad. The critical point of this plot is that the asymptotic corner of this summation curve occurs at 40 ms, which is the value of the time constant assumed for the neural signal in this example. Thus, the form of the BOLD amplitude summation series (Figure 4, column d) provides a direct empirical estimate of the time constant of neural integration down to the millisecond range. There is no limit in principle to the temporal resolution that can be achieved by this methodology since it is estimated from the amplitude variation of the BOLD signal as a function of stimulus duration, not from its temporal aspects.

The reduction in peak amplitude is captured in the red temporal summation curve of Figure 4Ad, showing a reduction by a factor of two for long durations. (The green curve in Figure 4Ad represents the values expected for fully proportional linear summation of the energy in the stimulus pulse; it is an accelerating curve due to the logarithmic abscissa.) The case of the full-rectification model of the metabolic demand is shown in Figure 4Ac. Here the second neural response peak (i.e., that from the stimulus offset) plays a key role in varying the BOLD waveform, which first extends in time and then shows a two-peaked structure with reduced amplitude for long-duration stimuli.

Figure 4Ba shows the half- and fully rectified version of the metabolic demand to the same pulse duration series, where the neural inhibition is now assumed to be reduced in energy by 1.5% relative to the excitation. This small imbalance is magnified by the convolution with the sustained stimulus, and thus it results in a sustained component that is 12% of the amplitude of the initial transient (blue curve in Figure 4Ba) and then into an almost fully sustained set of BOLD response functions (Figure 4Bb, compare to Figure 4Ab). Thus, the form of the BOLD response functions can be strongly diagnostic of even slight variations in the properties of brief neural signals. Moreover, the nature of the metabolic demand function (half- or fully rectified) has a big impact on the form of the BOLD response, determining whether or not an offset peak occurs at the tail of the responses even when they are sustained (Figure 4Bc). Such a peak has been reported in some studies but is not always evident. Thus it remains an empirical question to what extent rectification is representative of BOLD waveforms; intermediate forms of the rectification model are required to capture the empirical properties in detail.

Note that the amplitude series in Figures 4Bb,c show bands of denser packing of the functions, where the amplitude changes were not spaced in proportion to the doublings of stimulus duration. Viewed in terms of the sequence of BOLD waveforms in Figures 4Bb,c, the regions of dense packing form an intermediate “shelf” or partial asymptote in peak amplitude summation plots of Figure 4Bd. It is again evident that the onset of this intermediate shelf in the summation curve corresponds to the integration time of the underlying neural signal, while the second asymptote at higher amplitude corresponds to the ~5 s integration time of the MRK (HRF). Accurate measurement of such summation functions can therefore provide discriminative characteristics that, when interpreted through the non-linear model structure, can provide estimates of both the neural and the metabolic time constants in the neural-to-BOLD signal chain.

This point is emphasized by the response set in row c of Figure 4, which probes the effect of varying the time constant of the neural transient. The key difference from the parameters used in row b of Figure 4, is that the neural time constant was doubled from 50 to 100 ms (and the excitation/inhibition imbalance was also increased to 7% to maintain the same form of offset peak). It is evident that (i) the summation curve (Figure 4Cd) takes a measurably different form, and that (ii) the accuracy of estimation of the neural time constant is limited not by the BOLD time constant but only by the variability of the BOLD amplitude measures. For example, this analysis shows that the neural time constant is estimable to within about 0.1 log units if the BOLD response functions can be measured to an achievable accuracy of about 10%.

The final case of duration summation analyses (Figure 4Dd) shows the NDFO predictions of increasing the excitation/inhibition imbalance of the neural response to 50%, illustrative of a system that is predominantly sustained in nature. Under these conditions, the impact of the initial transient becomes essentially negligible, and the summation curves (Figure 4Dd) become indistinguishable from proportional summation (i.e., they run parallel to the green curve). This manipulation illustrates that the power of the NDFO analysis depends on the neural processing being predominantly transient, and that the properties of the underlying neural mechanisms would not be accessible to this form of analysis in predominantly sustained systems. Luckily, however, the well-established deviation from proportionality for short-duration stimuli implies that the neural system is, in practice, predominantly transient and is therefore amenable to this form of NDFO analysis.

To validate our conceptual breakthrough for the estimation of neural response dynamics from BOLD response measurements we performed a computational analysis of the implications of the BOLD fMRI data collected on a typical human subject with corrected to normal vision, a temporal modulation sensitivity function within normal limits, typical event-related potentials to contrast reversal (onset latency of 52 ms and N1 latency of 115 ms) and a typical BOLD response time course as evident from Figure 5 and a number of previous retinotopic and localizer scans. The analysis was conducted in a full-range temporal summation paradigm (Figure 5A) and provides not only the estimated neural response waveform but confidence limits on the waveform parameters that are sufficiently narrow to indicate that the results are within the expected range for early cortical response dynamics.

The responses from primary visual cortex (~250 voxels) for a 1-h session from a typical subject are shown in Figure 5B. They show a progressive growth in amplitude with stimulus duration from the shortest to the longest, but note that the amplitudes do not increase by a factor of two (the linear GLM prediction; cyan dashed curves in Figure 5B) as duration is increased. We therefore need a non-linear model to understand the processing mechanism for these effects.

The data have thus been fit through our NDFO procedure of a linear convolution of the stimuli from durations from 16.67 ms to 1.067 s in factors of 2 (Figure 5A) with each neural component kernel to generate the neural signals at each duration (Figure 5C), a non-linear coupling to generate the component metabolic demands, a linear summation over the component demands to generate the net metabolic demand, and a linear convolution with a MRK (Figure 5D) to generate the ultimate BOLD signal prediction (Figure 5B). Figure 5E is a comparison of the amplitudes of the fitted model waveforms with those of the data, as analyzed in Section “Discussion.” Optimization of the BOLD prediction (blue curves in Figure 5B) to the data determines the temporal features of the requisite neural signals making up each waveform, allowing us to specify the predominant time constants for the underlying neural signals. These features are defined by the following model population parameters for the net neural signal throughout each studied cortical region: transmission delay, gamma exponent, recovery time constant, sustained/transient ratio, sustained decay, and off/on ratio (six waveform parameters). The amplitude, transient integration time, and sustained/transient ratio are specified to worthwhile precision (see Appendix). The off/on ratio parameter controlling the relative amplitude of the off-responses is not well-constrained for the duration-range used in this data set (Figure 5), but we note that this parameter could also be brought under tighter control by adding stimuli of longer duration (e.g., 8 s). This approach should readily provide tight constraints on the off/on ratio because the BOLD activity corresponding to the off-responses is pushed to a time period that is essentially independent of that for the on response. This kind of approach to the off transient had been taken successfully by Fox et al. (2005) and Lingnau et al. (2009), for example.