NeuroImage 94 (2014) 203–215

Contents lists available at ScienceDirect

NeuroImage journal homepage: www.elsevier.com/locate/ynimg

Deconvolution of neural dynamics from fMRI data using a spatiotemporal hemodynamic response function K.M. Aquino a,b,c,⁎, P.A. Robinson a,c,d,e, M.M. Schira f,g,h,j, M. Breakspear b,h,i,j a

School of Physics, University of Sydney, New South Wales 2006, Australia Queensland Institute of Medical Research, Herston, Queensland 4006, Australia Brain Dynamics Center, Sydney Medical School Western, University of Sydney, New South Wales 2145, Australia d Cooperative Research Center for Alertness, Safety & Productivity, Notting Hill, Victoria 3168, Australia e Center for Integrated Research & Understanding of Sleep, Glebe, New South Wales 2037, Australia f School of Psychology, University of Wollongong, New South Wales 2522, Australia g Neuroscience Research Australia, Randwick, New South Wales 2031, Australia h School of Psychiatry, University of New South Wales, Sydney, New South Wales 2031, Australia i The Royal Brisbane and Women's Hospital, Brisbane, Queensland 4029, Australia j The Black Dog Institute, Sydney, New South Wales 2031, Australia b c

a r t i c l e

i n f o

Article history: Accepted 3 March 2014 Available online 12 March 2014 Keywords: fMRI Deconvolution Spatiotemporal BOLD HRF stHRF

a b s t r a c t Functional magnetic resonance imaging (fMRI) is a powerful and broadly used means of non-invasively mapping human brain activity. However fMRI is an indirect measure that rests upon a mapping from neuronal activity to the blood oxygen level dependent (BOLD) signal via hemodynamic effects. The quality of estimated neuronal activity hinges on the validity of the hemodynamic model employed. Recent work has demonstrated that the hemodynamic response has non-separable spatiotemporal dynamics, a key property that is not implemented in existing fMRI analysis frameworks. Here both simulated and empirical data are used to demonstrate that using a physiologically based model of the spatiotemporal hemodynamic response function (stHRF) results in a quantitative improvement of the estimated neuronal response relative to unphysical space–time separable forms. To achieve this, an integrated spatial and temporal deconvolution is established using a recently developed stHRF. Simulated data allows the variation of key parameters such as noise and the spatial complexity of the neuronal drive, while knowing the neuronal input. The results demonstrate that the use of a spatiotemporally integrated HRF can avoid “ghost” neuronal responses that can otherwise be falsely inferred. Applying the spatiotemporal deconvolution to high resolution fMRI data allows the recovery of neuronal responses that are consistent with independent electrophysiological measures. © 2014 Elsevier Inc. All rights reserved.

Introduction Functional MRI has progressed beyond “first-order” statistical inferences (classic task-locked activations) to uncover interactions between cortical regions (Friston et al., 2003; Horwitz, 2003). Such inferences require a solution to the inverse problem — i.e., estimation of neural activity from the BOLD signal. Several solutions to the inverse problem have been proposed. Notable examples include the direct deconvolution of fMRI data using a measured HRF (Gitelman et al., 2003; Glover, 1999) or the inversion of a forward model, such as implemented in dynamic causal modeling (DCM) (Friston, 2009). These existing analyses can then provide “second-order” insights into brain function, such as inferring interactions amongst cortical regions. ⁎ Corresponding author at: School of Physics, University of Sydney, New South Wales 2006, Australia. E-mail address: [email protected] (K.M. Aquino).

http://dx.doi.org/10.1016/j.neuroimage.2014.03.001 1053-8119/© 2014 Elsevier Inc. All rights reserved.

Despite the successes of these approaches (David et al., 2008; Friston, 2009; Stephan and Roebroeck, 2012), existing solutions are nonetheless limited as they typically treat the temporal dynamics of the hemodynamic response independently from its spatial dynamics. Such approaches miss intrinsic spatiotemporal hemodynamics (Gardner, 2010; Kriegeskorte et al., 2010) such as those that arise from draining veins (Bianciardi et al., 2010; Glover, 1999) and hemodynamic interactions between different vascular layers (Tian et al., 2010). Therefore, hemodynamic responses that couple space and time, such as hemodynamic waves (Aquino et al., 2012, 2014; Drysdale et al., 2010) cannot be adequately treated. Neglecting these spatiotemporal hemodynamics has important consequences for inferences of neural activity. For example time-lag based measures of effective connectivity, such as Granger causality, are vulnerable to spatiotemporal variability in the HRF (Chang et al., 2008). The basis of these effects on inferences of cortical connectivity is that spatially distributed delays may be misattributed to neural dynamics. Therefore such techniques rely upon an accurate

K.M. Aquino et al. / NeuroImage 94 (2014) 203–215

spatiotemporal characterization of the hemodynamic effects that appear in the BOLD signal. Recent research has shown that the spatiotemporal hemodynamic response function (stHRF) contains complex spatiotemporal structure (Aquino et al., 2012; Bießmann et al., 2012). To date, fMRI analysis procedures have explicitly or implicitly employed an HRF where space and time are described independently, typically assuming that the spatial profile to be a Gaussian function in three dimensions; the temporal profile is typically described by a double gamma function [e.g. (Boynton et al., 1996)] or using output from the balloon model (Buxton et al., 2004). We argue that fMRI analysis could be significantly improved with an stHRF that is spatiotemporally non-separable and derived from physiology. We note three key aspects of recent studies: (i) Inferences of specificity and localization of brain activation detected in fMRI improve when analyzing activation on the brain's “natural” coordinates embedded on the cortical surface (Chen et al., 2011), as opposed to analyzing in naive voxel space. (ii) According to a recent model and corresponding experiments, spatiotemporal hemodynamics exhibit traveling waves (Aquino et al., 2012, 2014) thus motivating the need to analyze fMRI using tools from wave theory. (iii) Advancing the utility of fMRI will come from signal analysis techniques that couple space and time (Kriegeskorte et al., 2010). A recent data driven study (Bießmann et al., 2012) showed, that using a proposed nonseparable stHRF improves localization relative to a separable approximation. The purpose of the present study is to estimate spatiotemporal neural activity by deconvolving functional magnetic resonance imaging (fMRI) data using an stHRF derived from physiology. In particular, this study aims to refine estimates of neural activity by solving the fMRI inverse problem using an stHRF – on the cortical surface – that is based on a poroelastic treatment of cortical tissue (Aquino et al. 2014; Drysdale et al., 2010). We first show how spatiotemporal fMRI data can be inverted by using a deconvolution procedure that relies on the Wiener filter with an stHRF. We then compare predictions of the BOLD response with the model-derived stHRF against the previously assumed space– time separable HRF. We then test our deconvolution methods on simulated data. Finally, we use these results to deconvolve fMRI data we recently acquired in an evoked study of the visual cortex (Aquino et al., 2012) and hence show that using an stHRF based on physiology improves estimates of neural activity relative to the unphysical assumption of space–time separability (it implies infinite signal propagation). These improvements include more accurate estimates of the spatial width of neural activation and avoidance of confounding estimates of neural activity. Methods Deconvolution of fMRI responses Changes in neural activity drive a series of physiological processes that lead to changes in the BOLD signal. Although the details of these processes are not completely understood, their effect can be summarized via the combined hemodynamic response to a discrete input — the HRF. As long as the BOLD response to a neural drive is linear, then the resultant BOLD response Y can be treated as a convolution of the HRF G and neural drive ζ Y ðr; t Þ ¼ Gðr; t Þ⊗ζ ðr; t Þ:

Y ðk; ωÞ ¼ Gðk; ωÞζ ðk; ωÞ;

ð3Þ

where the arguments denote wavenumber k and frequency ω. Therefore, any convolution can be written as a product of the source ζ and the HRF G in Fourier space via the Fourier convolution theorem (Melrose and McPhedran, 1991). We can deconvolve Y to obtain an estimate of the neural drive ζ^ in Fourier space by calculating the ratio: Y ðk; ωÞ ζ^ ðk; ωÞ ¼ ; Gðk; ωÞ

ð4Þ

and then by taking the inverse Fourier transform to obtain an estimate ζ^ in coordinate space: ζ^ ðr; t Þ ¼

3

e ∫ ðd2πkÞ ∫ dω 2π

−iðωt−krÞ

3

Y ðk; ωÞ ; Gðk; ωÞ

ð5Þ

where the integrals are over all k and ω. Weiner filtering When implementing Eq. (5) an issue that arises is the presence of singularities (or near-singularities) arising from zero (or small) values in the denominator G. These zero or small values in G are evident at high spatial and/or temporal frequencies. This is because G is analytically described or measured in a controlled way i.e., due to specific the evoked response as opposed to the BOLD response to a series of stimuli, hence it has low noise levels with the effect that G(k, ω) rapidly decays to zero with increasing k or ω. However, since Y typically contains physiological [e.g., Frank et al. (2001)] or scanner noise (Krueger and Glover, 2001), Y(k, ω) typically decays to zero more slowly than G. Hence the ratio (5) becomes large, in essence amplifying noise in the higher frequencies. This is demonstrated in Fig. 1, for a one dimensional response (described in Deconvolution of synthetic data), by comparing the power spectrum, Z

2

jGðωÞj ¼

þ 

2

jGðk; ωÞj dk;

ð6Þ

with the corresponding |Y(ω)|2. Above ≈ 0.1 Hz, the contribution due to G(ω) decays more steeply than |Y(ω)|2 so that the estimate (5) for ζ^ amplifies high frequency components. Since the hemodynamic response acts as a low-pass filter with a cutoff of 0.1 Hz (Robinson et al., 2006; Zarahn et al., 1997) frequencies above ≈ 0.1 Hz are considered to be dominated by artifacts, including: aliasing from movements due to breathing (Frank et al., 2001), aliasing from cardiac pulsation (Frank

|G(ω )|2 |Y (ω )|2 |Y (ω )|2/ |G(ω )|2 10−5

ð1Þ

This convolution, represented as ⨂, is given by

 ′ ′ ′ ′ Gðr; t Þ⊗ζ ðr; t Þ ¼ ∫d r ∫dt G r−r ; t−t ζ r ; t : 3 ′

space and in time. By taking the Fourier transform we can rewrite this convolution as

Power

204

′

10−10 0

ð2Þ

where the integrals are over all space r′ and time t′. This framework treats the BOLD signal as a filtered version of the neural drive in both

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

f (Hz) [= ω/(2π)] Fig. 1. Amplification of high frequency noise by naive deconvolution. The figure compares the average power spectra of the HRF |G(ω)|2, the BOLD signal |Y(ω)|2, and the ratio of these spectra |Y(ω)|2/|G(ω)|2. The functions G and Y are calculated in stHRF derived from physiology.

K.M. Aquino et al. / NeuroImage 94 (2014) 203–215

et al., 2000; Wang et al., 2005), and thermal measurement noise in the fMRI scanner (Krueger and Glover, 2001). There are various schemes to account for inherent noise contributions. Of particular interest is the employment of a filter whose design is informed by the likely noise structure of the signal. For instance, we can modify Eq. (1) to include a noise term N(r,t): Y ðr; t Þ ¼ Gðr; t Þ⊗ζ ðr; t Þ þ Nðr; t Þ:

ð7Þ

^ similar to Eq. (4), as the product of Y and a We seek an estimate for ζ, factor D, ζ^ ðk; ωÞ ¼ Dðk; ωÞY ðk; ωÞ:

ð8Þ

A solution in the least squares sense is the Wiener filter Dðk; ωÞ ¼

1 jGðk; ωÞj2  ;  Gðk; ωÞ Gðk; ωÞj2 þ NSRðk; ω

ð9Þ

where NSR represents the noise to signal ratio 2

NSRðk; ωÞ ¼

jNðk; ωÞj : jSðk; ωÞj2

ð10Þ

The signal term S represents the output of the deconvolution. A simple ansatz is to assume that the noise is white; i.e., independent of time and space, and that the signal is a delta function in space and time hence yielding: 2

NSRðk; ωÞ ¼ σ ;

ð11Þ

where σ2 is a constant that can be estimated from measurements of the signal in the absence of experimental manipulations. In summary, the deconvolution procedure is as follows: (i) Determine the HRF as shown below; (ii) Take the Fourier transform of the BOLD response Y and the HRF G; (iii) Estimate the noise to signal ratio; and (iv) Use these functions to estimate the neuronal drive ζ^ . The spatiotemporal hemodynamic response function The neuronal estimate derived from deconvolution is dependent on the nature of the HRF employed. Here we detail the two different response functions, one derived from previous empirical estimates, and an stHRF derived from a physical theory of hemodynamics. Separable assumptions of the stHRF As reviewed above, the spatiotemporal hemodynamic response function has been generally assumed to derive from a product of a spatially varying component with a temporal hemodynamic response.

(a) 20

This time course corresponds to the evoked temporal BOLD response that has been well studied and is treated as either an empirical fit of the evoked response; e.g., by the sum of two gamma functions (Handwerker et al., 2004), or with physiological models of hemodynamics, namely the Balloon model (Buxton et al., 1998, 2004; Friston et al., 2000). The spatial component is often approximated by a Gaussian spatial kernel on the basis that: (i) the spatial point spread functions in BOLD are approximately Gaussian (Engel et al., 1997; Menon and Goodyear, 1999; Parkes et al., 2005; Shmuel et al., 2007), (ii) this is the simplest statistical ansatz (Worsley and Friston, 1995); and/or (iii) the dispersion of vasodilatory gases, such as nitric oxide, are likely to be spatially diffusive (Friston, 1995). Together, these properties correspond to the following mathematical form for a separable stHRF denoted by Gs,
 2 2 GS ðr; t Þ ¼ exp −jrj =Δr HRFðt Þ;

ð12Þ

where Δr denotes the spatial spread, the time course HRF(t) is the temporal hemodynamic response, r are the spatial three dimensional coordinates and t is time. A consequence of having this mathematical form for G is that if the BOLD responses shows spatially distributed delays, these will be attributed to neuronal delays, as this hemodynamic response function does not include terms that relate r and t. Here the HRF(t) is calculated with the balloon model described in (Friston et al., 2000), using the parameters below. A one dimensional example is visualized in Fig. 2a. We note that this space time separability approximation leads to an unphysical description of the hemodynamic response. The stHRF Gs theoretically implies signal propagation, that leads to spreading of the HRF, occurs instantaneously and practically it implies that the signal propagation is much faster than the sampling rate of fMRI scanners. However, this spreading is determined by the hemodynamic transit time (≈ 1 s) over the ≈ 1 mm spacing between arteries and veins, thus resulting in speeds of order 1 mm s−1 (Aquino et al., 2012, 2014). stHRF derived from physiology We now consider the stHRF derived from a spatiotemporal model of hemodynamics, treating cortical tissue as a porous medium of dense vasculature that penetrates the cortical tissue in different directions. This is a mean field approximation of cortical tissue and is used to model spatiotemporal changes in blood volume density, concentration of oxygenated blood, mean blood fluid velocity and pore pressure. This theory models the changes in cortical tissue due to increases in blood flow constrained by physical laws — notably conservation of mass, momentum, and total hemoglobin. These laws allow the model to be formalized as a closed set of nonlinear partial differential equations that are driven by neural activity.

(b)

GP

205

20

GS

0.8

0.8 15

0.6 0.4

10

0.6

t (s)

t (s)

15

10

0.4

0.2

0.2

5

5 0 0 −5

0

x (mm)

5

−0.2

0 0 −5

0

5

x (mm)

Fig. 2. Spatiotemporal hemodynamic response functions. (a) Space–time separable ansatz Gs due to a localized neural drive using σ = 3 mm. (b) Spatiotemporal HRF derived from physiology Gp. The full parameters for these stHRFs are shown in Appendix A. The color scale indicates the BOLD amplitude, normalized to have maximum amplitude of unity hence the responses are dimensionless.

206

K.M. Aquino et al. / NeuroImage 94 (2014) 203–215

By linearizing the dynamical equations and averaging through the thickness of the cortex (Aquino et al. 2014), the stHRF is calculated as a function of position on the cortical surface. Through this linearization, the BOLD response can be represented in Fourier space as the transfer function that links the BOLD response Y to a neural drive ζ:
TYζ ðk; ωÞ ¼

 −iω þ D=ρ f ðk2 −k3 Þ

k2 ν 2β þ k2z ν2β −ω2 −2iΓω 1  −ðω þ i12κ Þ2 þ ω2f h i9 8 −1   < k1 þ k2 −V 0 iω−C Z τ ðβ−2Þ þ η = :  1− : ; k2 −k3 −iω þ η þ τ−1

ð13Þ

where the full parameters are detailed in Appendix A. The transfer function describes the change in Y per unit change in ζ at the same k and ω. Hence the linear BOLD response to a neural drive ζ can be calculated in Fourier space as, Y ðk; ωÞ ¼ T Yζ ðk; ωÞζ ðk; ωÞ:

ð14Þ

This can be calculated in co-ordinate space as

Y ðr; t Þ ¼

∫

d2 k ð2π Þ2

∫

dω −iðωt−krÞ T ðk; ωÞζ ðk; ωÞe : 2π Yζ

ð15Þ

The stHRF Gp is calculated as the BOLD response to a discrete and impulse neural drive, i.e., a delta function δ(r,t), hence substituting into Eq. (15) yields, 2

GP ðr; t Þ ¼

T ∫ ðd2πkÞ ∫ dω 2π

−iðωt−krÞ : Yζ ðk; ωÞe

2

∫ ðdk ∫ dω T 2πÞ 2π x

−iðωt−kx xÞ : Yζ ðkx ; 0; ωÞe

value of these parameters was estimated from previous studies and are shown in Appendix A, and the same parameters were used for both the stHRFs. The parameters unique to each stHRF are as follows: for Gs, Δr, and for Gp, (kz, ρf, CZ, νβ, Γ, D). We note that the parameters kz,Cz are functions of the cortical thickness L, and ρf is the mass density of blood, hence these parameters are constant. In addition, we note that D and Γ are related, thus together we reduce the variable parameter set of Gp to (Γ,vβ). Each of these parameters can be inferred from experiment or determined independently as they are based on physiological parameters (Aquino et al., 2012, 2014). Generation of synthetic BOLD data The deconvolution procedure will be first tested on synthetic BOLD data generated by the convolution of a neural drive, ζ, with the hemodynamic response function G. The neural drive is modeled as a simple function in time and in space, particularly in one dimension x as ζ ðx; t Þ ¼ e

ζ ðx; t Þ ¼

ð17Þ

Visualization of the one dimensional stHRF [Eq. 17] in Fig. 2b, shows the traveling wave character of this stHRF. The propagation of the BOLD response is damped as it travels left and right of the stimulus located at x = 0, t = 0. Physically, this BOLD response is generated by key stages in the blood volume dynamics: (i) The increase of the neural drive induces vasodilation, causing an increased influx of blood increasing local vessel volume and pressure. (ii) Blood flows to regions of lower pressure, redistributing pressure throughout the medium. These pressure changes induce blood volume and fluid velocity in adjacent cortical tissue that propagate at a speed vβ. (iii) As these coupled changes propagate, their influence is dampened by processes that include viscosity and drainage of blood. This results in damped hemodynamic waves, and is parametrized by Γ. In addition, changes in local deoxyhemoglobin (dHb) content occur. Lowering of dHb concentration is caused by the influx of oxygenated hemoglobin (oHb) by arterial blood, and the removal of dHb by local outflow. During this process the concentration of dHb is increased by the passive diffusion of oxygen into the cortical tissue, converting oHb to dHb. Together, these processes result in an initial decrease of local dHb concentration followed by a delayed increase. The combined changes of blood volume and dHb lead to change in the BOLD signal.

−x2 =σ 2x

h i 2 2 exp −ðt−t 0 Þ =σ t ;

ð18Þ

where σx is the spatial spread chosen for ζ to have full width half maximum of 0.5 mm, σt is the temporal spread chosen for ζ to have full temporal width half maximum of 0.5 s, and t0 is the onset time set at 2 s. Using Eq. (3) we can generate BOLD responses using the two different stHRFs. Additionally, we model the neural drive of two neural sources separated by distance d and centered at ±d/2 as

ð16Þ

Equivalently, the one dimensional stHRF is calculated as,

GP ðx; t Þ ¼

Parameter estimation of the stHRFs The combined parameter set of the two stHRFs contains overlapping parameters that allows direct comparison between the two stHRFs. The overlapping parameters are comprised of the following h
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi i set: k1 ; k2 ; k3 ; α ð¼ 1=βÞ; κ; γ; ω f ¼ γ2 −κ 2 =4 ; V 0 ; τ; E0 ð¼ τηÞ . The

h i 2 2 2 2 1 2 2 −ðxþd=2Þ =σ x −ðx−d=2Þ =σ x exp −ðt−t 0 Þ =σ t e þe : 2

ð19Þ

Functional MRI data The deconvolution will be then employed to analyze fMRI data collected from the visual cortex, following presentation of discrete spatial stimuli (Aquino et al., 2012). The data were acquired in a 3 T fMRI Phillips Achieva scanner using a 8 channel head coil to resolve fMRI with 1.5 mm cubic voxels and a TR of 2 s with coverage restricted to the occipital pole. The visual stimulus was an image of three rings at 0.6°, 1.6°, and 3° eccentricity that oscillated in contrast at a rate of 2 Hz for 8 s during the “on” phase, and was not present during the “off” phase that lasted 12 s. These rings were presented simultaneously, and were retinotopically far enough such that the responses from the first ring (which we analyzed) did not interfere with the outer two rings, as described in Aquino et al. (2012). These rings induced approximately straight lines on the visual cortex, which we exploit in our analysis (see below). Furthermore, the presentation of this stimulus was temporally jittered to yield an effective temporal sampling of the response at 250 ms resolution. The experiment was repeated on four healthy subjects. The BOLD responses were projected onto the cortical surface using a surface reconstruction from a gray and white matter segmentation of a T1 weighted image sampled with 0.75 mm isotropic voxels. Because the BOLD response to the stimulus induced approximately straight lines on the cortical surface, the hemodynamic response depends only on time t and the perpendicular distance x from that neural response locus. We estimated the location of the centerline of the primary response to be at 0.6° on a flattened representation of the cortical surface and then measured the average BOLD signal at various distances orthogonal to

K.M. Aquino et al. / NeuroImage 94 (2014) 203–215

the synthetic BOLD response appears as a simple low pass filter of space and time with a temporal profile independent of distance. The simulations demonstrate that although the spatial extent of the expected BOLD response might be similar under these two stHRFs, their spatiotemporal dynamics are markedly different. In addition, the existence of wave dynamics in Gp means that the hemodynamic responses from nearby sources of neural activity may show wave–wave interactions. In Fig. 4 we use Eq. (19) to calculate the BOLD response due to two neural drives whose centers are separated by 6 mm, as shown in Fig. 4(a). While the response of these sources calculated using Gs is the sum of two Gaussian functions, as in Fig. 4(c) and (d), when Gp is used the BOLD responses can interact to generate an additional response between the two sources, as shown in Fig. 4(b). The differences between these two responses are highlighted in the time to peak profile in Fig. 4(e). The BOLD profile derived using Gs is flat — indicating the same temporal response throughout, whereas the profile derived using Gp suggests that there is a local response that peaks at t ≈ 6 s at x = ± 3 mm, and an additional BOLD response at x = 0 mm that peaks at 7 s.

this response, restricted to V1, and as a function of time since stimulus onset. As a result each hemodynamic response, Y(x,t) was sampled across 10 mm, along the interval x = [−5,5] mm in steps 0.75 mm, and across 20 s with time axis t = [0,20] s in steps of 0.25 s.

Results The results in this section are fourfold: we first describe the properties of the simulated BOLD response due to a localized neural drive; we then demonstrate the hemodynamic interactions due to two neural drives; we next detail the results from the deconvolution of these synthetic sources; and finally, these outcomes guide deconvolution of experimental fMRI data in visual cortex.

Simulated BOLD responses: the effect of stHRFs Given a neural drive, the properties of the simulated BOLD responses depend on the characteristics of the stHRF. In Fig. 3(a) we show a short and relatively local neural drive ζ. By using the physiological stHRF Gp, we show hemodynamic traveling waves traveling outward from x = 0 as seen in Fig. 3(b), with speed vβ and damping at a rate Γ. The spatial extent of this response is captured by measuring the maximum amplitude of the BOLD response as a function of distance, as seen in Fig. 3(c), whereby the BOLD response covers a wider range than the neural drive. The properties of the hemodynamic waves are presented in measures of the time to peak as a function of distance, for example Fig. 3(d) reveals time delays left and right from x = 0 due to activity traveling away from the central region of the neural drive. In contrast, the BOLD response using the separable stHRF ansatz Gs [Fig. 3(b)] does not display the temporal delays as a function of distance, as evidenced in Fig. 3(d). Since Gs is separable in space in time,

ζ

(a)

(b)

Deconvolution of synthetic data Here, we validate the deconvolution procedure on synthetic data generated from the spatiotemporal hemodynamic model. The parameters (vβ, Γ) for Gp were taken from the mean parameter estimates from our companion study (Aquino et al., 2012) and the parameter Δr for Gs was guided by typical measurements from the point spread function in 3 T (Parkes et al., 2005). The parameters common to both responses were taken from independent studies and are provided in Appendix A. Testing a simple neural drive We deconvolve the BOLD signal calculated from a single neural drive shown in Fig. 3(b). This requires an estimation of the NSR. Since we have

(c)

GP

20

20

5

0.2

10

0.4 0.2

5

t (s)

t (s)

15

0.6

0.6 0.4

0.8

0.8 15

10

GS

20

0.8 15 t (s)

207

0.6 0.4

10

0.2 5

0

0

0 −5

0 x (mm)

5

0 −5

0 x (mm)

0 −5

−0.2 0 x (mm)

5

(e)

(d) 1

10

0.8

8

time to peak (s)

max amplitude

5

0.6 0.4

4 2

0.2 0 −5

GP GS ζ

6

0

x (mm)

5

0 −5

0

5

x (mm)

Fig. 3. Dependence of the BOLD response on the stHRF. (a) Localized neural drive ζ calculated using Eq. (18). (b) BOLD response calculated with Gp. (c) The BOLD response calculated with Gs. (d) Maximum amplitude of the responses as a function of distance. Red line is ζ, the black is the BOLD response using Gp and the blue is the BOLD response using Gs. (e) Time for the response to peak vs. distance, using the same convention as in (c). Each response has been normalized to have a maximum amplitude of unity, and hence dimensionless.

208

K.M. Aquino et al. / NeuroImage 94 (2014) 203–215

ζ

(a)

(b)

(c)

GP

20

20 0.8

10

0.4 0.2

5

0.2

t (s)

t (s)

t (s)

5

15

0.6

0.6 0.4

0.8

0.8 15

15 10

GS

20

0.6 0.4

10

0.2 5

0

0

0 −5

0 x (mm)

0 −5

5

0 x (mm)

0 −5

−0.2 0 x (mm)

5

(e)

(d) 1

8

time to peak (s)

0.8

max amplitude

5

0.6 0.4 0.2 0 −5

0

6

4

2

0 −5

5

x (mm)

GP GS ζ

0

5

x (mm)

Fig. 4. Dependence of the BOLD response on the stHRF due to two neural drives. (a) The neural drive ζ with two drives separated by 6 mm and centered at x = ±3 mm using Eq. (19). (b) BOLD response calculated with Gp. (c) BOLD response calculated with Gs. (d) Maximum amplitude of the responses as a function of distance, where red line is ζ, the black is the BOLD response using Gp and the blue is the BOLD response using Gs. (e) Time for the responses to peak as a function of distance, using the same convention as in (c). Each response has been normalized to have a maximum amplitude of unity, and hence dimensionless.

the analytic form of the neural drive, we can quantify the difference between ζ^ to the original drive ζ with the following difference metric ϵ, h i2 ∫dx∫dt ζ ðx; t Þ−ζ^ ðx; t; NSRÞ h ϵðNSRÞ ¼ 2 2 ∫dx∫dt ζ ðx; t Þ þ ζ^ ðx; t; NSRÞ ;

ð20Þ

where ζ^ ðx; t; NSRÞ is the deconvolved estimate of the neural drive for a given NSR. This metric varies from 0 to 1, where the highest values represent the greatest difference. We use this metric to deconvolve the BOLD signal using both Gp and Gs. The different results achieved from these are shown in Fig. 5. We find that NSR ≈ 200 minimizes ϵ, for estimates using Gp or Gs, hence optimizes the deconvolution procedure. Note that ϵ plateaus for NSR values too high (≳ 2 × 104) or too low (b 1 × 10−3), and that for Gp, ϵ is lower relative to the estimate from Gs implying that the signal is best deconvolved with the signal from which it is generated, namely Gp.

G

P

1

G

S

ε(NSR)

0.8 0.6 0.4 0.2 0 10−10

10−5

100

105

1010

NSR Fig. 5. The difference metric ϵ as a function of NSR for neural drive estimates ζ^ with Gp (in solid) and Gs (in dashed lines).

In practice however, the NSR needs to be determined by a data drive heuristic. To achieve this, we note that the Wiener filter has a low-pass character with a cutoff frequencies (fc,kc) that are roughly determined the highest frequencies (k, ω) that satisfy |G(k, ω)|2 N NSR. In the temporal domain, due to the low-pass character of the BOLD signal relative to neural activity, fc ≈ 0.1 Hz. Spatially, we use the cutoff at kc/(2π) ≈ 480 m−1 corresponding to wavelengths of ≈2 mm. We first demonstrate the deconvolution estimates using Gp. Exemplar deconvolutions with a range of frequency cutoffs are shown in Fig. 6. The various deconvolution estimates seen in Figs. 6(b)–(d) recover the broad properties of the original neural drive, as seen in Fig. 3(a), particularly that ζ^ is Gaussian in space and time. We have thus successfully deconvolved the traveling wave behavior seen in the BOLD response. Differences in the deconvolved estimates depend on changes in frequency cutoffs. For instance, a high frequency cutoff (Fig. 6b) leads to a signal that has a similar spatial and temporal extent of the neural drive relative to estimates with a lower cutoff. This comes at the cost of high frequency noise in the estimate. Lowering this cutoff in fc and kc leads to a smoother signal at the cost of over estimating the spatial extent and the temporal duration [see Figs. 6(c) and (d)]. These features are shown in cross-sections in time and space in Fig. 7. In addition, the cross sections in Figs. 7(a) and (b) show that the estimates for ζ^ contain negative components – after t = 5 s, and beyond |x| N 2.5 mm – which are not present in ζ. These oscillatory artifacts are due to the Wiener filtering: Higher frequencies are filtered such that the estimate ζ^ is dominated by larger oscillations that become more pronounced with lower frequency cutoffs as seen in Fig. 7. As discussed above, typical BOLD signals contain experimental artifacts arising from scanner noise and endogenous physiological activity. ^ we repeat the To explore how these effects influence the estimates of ζ, above analysis after adding white noise at a level of 5% of the peak simulated BOLD amplitude. The deconvolved responses are shown for two cut-off frequencies in Fig. 8. With lower frequencies the estimate for ζ^ is smoothed for lower cutoff frequencies (Fig. 8b), and in the processes

K.M. Aquino et al. / NeuroImage 94 (2014) 203–215

(a)

(b)

20

0.5

209

(c)

20

(d)

20

0.8

20 0.8

0

15

15

15

0.6

10

0.4

−0.5

10

10

5

−1

5

0 −5

−1.5

0 −5

0 x (mm)

0.2

t (s)

t (s)

t (s)

0.4

0 x (mm)

0.6

10

0.4 0.2

5

5 0

0 0 −5

5

15

0.2

0 −0.2

5

t (s)

0.8 0.6

0 x (mm)

0 −5

5

0 x (mm)

5

Fig. 6. Estimates of the neural drive using deconvolution. Panels (a)–(d) are deconvolved estimates of ζ^ using (a) fc = 2 Hz, kc/(2π) = 1000 m−1 (b) fc = 0.1 Hz, kc/(2π) = 500 m−1, (c) fc = 0.1 Hz, kc/(2π) = 200 m−1 and (d) fc = 0.05 Hz, kc/(2π) = 200 m−1. Each response has been normalized to have a maximum amplitude of unity, and hence dimensionless.

overestimates the spatial and temporal extent. The general properties of ζ^ are still recovered under the range of frequency cutoffs we previously specified, showing that the cutoffs fc = 0.1 Hz, kc/(2π) = 500 m -1 are sufficient to ensure that the spatiotemporal extent of the neural drive is obtained. In standard deconvolution techniques, inverting signals using the same function that generates the signal optimizes the resulting inversion estimate (Gitelman et al., 2003; Kerr et al., 2009; Lee et al., 1995). To demonstrate the artifacts that a “mismatched” stHRF generates, we deconvolved our synthetic data with the separable stHRF Gs shown in Figs. 9(a)–(d). These estimates on average do not recover ζ as well as Gp as evidenced by the larger value of the difference metric ϵ in Fig. 5. Furthermore, we note three important artifacts: (i) The spatial and temporal extent of ζ^ are more sensitive to frequency cutoffs compared with the estimates using Gp. (ii) For a range of frequency cut-offs ζ^ contains nonexistent neural sources. For example in Fig. 9(c) there is a strong feature that occurs from 5 to 10 s at x ± 2.5 mm, and (iii) For low frequency cut-offs i.e., Fig. 9(d), ζ^ contains waves as indicated by the time to peak as a function of distance (in dashed-black) indicating propagation of the signal left and right from x = 0 which is not present in ζ. For completeness, if we repeat this analysis with Gp as the mismatched filter (i.e. BOLD data is generated by Gs), the artifacts are less severe, and mainly result in the overestimation of the spatial width of the neural drive. These results are presented in Appendix B. An important issue when using any model is the balance between the goodness of fit of the model and its complexity. When deconvolving simulated data (where the ground truth neural drive is known), this

balance can be formally estimated by comparing the ζ^ estimates, generated with Gp and Gs, to the original drive ζ using the Bayesian information criterion BIC (Schwarz, 1978) - this rewards models by their variance explained and penalizes model complexity. The best model for a given data set has the lowest BIC value. Here, we calculated the BIC in Fig. 10 with 13 and 10 parameters for Gp and Gs respectively across a broad range of NSR. These analyses indicate that Gp outperforms Gs – despite its added complexity – when the data is generated by Gp for a wide range of NSR values. For extremely low NSR values, Gs is the preferred model because the inversion estimate is noisy (due to lack of filtering) and the variance that cannot be explained is high for both models. According to Occam's Razor, the least complex model Gs should be preferred in this noisy setting. These results further indicate that Gp will induce specific effects that need to be accounted for in deconvolution. Deconvolution of multiple neural sources We explore how the estimates of ζ^ with multiple neural drives are influenced following deconvolution. The neural drive was modeled by the sum of two Gaussians as in Eq. (19), with a separation of d = 6 mm as shown in Fig. 11(a). The predicted BOLD response is shown in Fig. 11(b). As expected, using Gp we recover ζ as shown in Fig. 11(c). This reconstruction contains a low amplitude artifact at x = 0 mm from 4 to 8 s due to the suppression of high frequencies, due to Weiner filtering, that would allow a rapid decay to zero. However, using the mismatched filter Gs we note that while we recover two neural sources centered at x = ± 4 mm as seen in Fig. 11(d), this inversion predicts a “ghost” neuronal response centered at x = 0 mm from ≈ 5–10 s. Deconvolution of empirical fMRI data

0.5

0.5

ζ (2,1000) (0.1,500) (0.1,200) (0.05,50)

We now apply the deconvolution procedure to the evoked spatiotemporal BOLD responses, studying the response in V1 to a discrete visual stimulus. The spatiotemporal BOLD response in an exemplar

(a)

ζ(0,t )

1

(b) 20

20

>

>

Average ζ

(b) 1

0

0.8

0.8 15

15

0

0.6

t (s)

0.6 10

0.4

t (s)

(a)

10

0.4 0.2

0.2

−0.5 −5

0

x (mm)

5

−0.5

5

5

0

10

t (s)

Fig. 7. Cross sections of the estimates of ζ^ described in Fig. 6. (a) Spatial profile of ζ^ averaged over the first 5 s. (b) Time variation on the centerline, x = 0. The different curves represent deconvolution estimates with frequency cutoffs fc and kc/(2π) in the pair [fc,kc/(2π)]. Each response has been normalized to have a maximum amplitude of unity, and hence dimensionless.

0

0

20 0 −5

0 x (mm)

5

−0.2

0 −5

0 x (mm)

5

−0.2

Fig. 8. The effects of noise in estimating ζ^ with two different frequency cutoffs. (a) Estimate for ζ^ using frequency cut offs fc = 0.1 Hz,kc/(2π) = 500 m−1 with Gp. (b) Estimate for ζ^ using frequency cut offs fc = 0.1 Hz, kc/(2π) = 200 m−1 with Gp. Each response has been normalized to have a maximum amplitude of unity, and hence dimensionless.

K.M. Aquino et al. / NeuroImage 94 (2014) 203–215

(a)

(b)

20

(c)

20

(d)

20

0.8

0

−2

5

10 5

−3

0 −5

0 x (mm)

5

0

−4

0 −5

15 t (s)

10

0.5

15

−1 t (s)

t (s)

15

0.2

10

0 x (mm)

5

0.6

15

0.4

0.4 10 0.2

0 5

−0.5

20

0.6

t (s)

210

0 −5

−0.2 −0.4 0 x (mm)

5

5 0 −5

0 0 x (mm)

5

−0.2

Fig. 9. Estimates of the neural drive by deconvolving with a mismatched filter Gs. Panels (a)–(d) are deconvolved estimates of ζ^ using (a) fc = 2 Hz, kc/(2π) = 1000 m−1 (b) fc = 0.1 Hz, k/(2π) = 500 m−1, (c) fc = 0.1 Hz, kc/(2π) = 200 m−1, and (d) fc = 0.05 Hz, kc/(2π) = 200 m−1. Each response has been normalized to have a maximum amplitude of unity, and hence dimensionless.

subject is shown in the Fig. 12(a). To deconvolve this dataset, we determine the stHRF by using the corresponding parameter estimates of velocity and damping obtained in our companion study (Aquino et al., 2012) with vβ = 2 mm s−1, Γ = 0.8 s−1. We use our previous frequency cutoffs used in testing for Gp. The result of deconvolving this response with Gp is shown in Fig. 12(b). The estimated neural drive ζ^ is localized around the centerline with a full width at half maximum (FWHM) of ≈ 5 mm, as seen in Fig. 13(a). This lasts for approximately the stimulus duration of 8 s. Because the hemodynamic model predicts hemodynamics that couple space and time, with traveling waves, the delays present in the BOLD response, shown in Fig. 13(b), are hemodynamic in origin. That is, Gp deconvolves such delays such that they are not apparent in ζ^ . We then deconvolved this same empirical BOLD data using Gs [Fig. 12(c)]. We find that although this estimate reveals a similar spatial [Fig. 13(a)], and temporal profile [Fig. 13(b) at x = 0], the combined spatiotemporal properties differ from the estimates obtained using Gp. The key difference is seen in the time to peak profile in Fig. 13(b) which implies traveling neural waves, propagating with velocity ≈ 1 mm s−1. This is at least 1–2 orders of magnitude slower than the corresponding neural waves measured independently (Benucci et al., 2007; Grinvald et al., 1994; Xu et al., 2007) implying that using Gs yields improbable spatiotemporal properties of the neural drive. The spatial extent estimated by both stHRFs are slightly greater than previously determined values (Aquino et al., 2012) and slightly larger than physiological estimates (Vanduffel et al., 2002). This discrepancy can be attributed to the smoothing artifacts due to the Wiener filter that induce enhancement of low spatial and temporal frequencies. There is a ringing artifact from 18 to 20 s due to Wiener filter. This particular artifact is due to the wrap around effects imposed by using the Fourier transform, and the filtering of high frequencies. Since the evoked response is zero-padded from [− 20,0] seconds in the deconvolution

−1

x 105

BIC

−1.5

−2

−2.5

G

P

G

S

−3 10−5

100

105

NSR Fig. 10. Model comparisons, calculated with BIC, for the estimate of the neural drive compared to the original drive as a function of NSR. BOLD data were simulated using Gp and deconvolved using Gp (solid line) and Gs (dashed line).

procedure, this imposes that the neural estimate needs to be a signal that will decay the BOLD response to zero at the boundary at t = 20 s. Since high frequencies are filtered out, this neural signal cannot decay sharply and hence causes an oscillation clearly noticeable from 18 to 20 s. This range was excluded from our time to peak analyses such as in Fig. 13. The deconvolution procedure was repeated for the remaining subjects, and is shown in the Supplementary text. The mean peak amplitude and mean time to peak as a function of distance are visualized in Fig. 14. We fitted the time to peak profiles with polynomial fits of order that minimized the Akaike Information Criterion (Akaike, 1974). From these fits we see that on average the traveling waves were deconvolved in Gp as noted by relatively uniform time to peak profile in Fig. 14(b). In contrast in the average the estimates from Gs are fitted best with a quadratic fit about x = 0 indicative of traveling waves [c.f. Fig. 3(e)]. We observe that time to peak profile is appears noisier with Gp. We determine that it is due to the fact that Gp deconvolves the waves away, leaving the localized neural response at the center. The remaining component of the neural drive is the residue that contains the peripheral response, which is thus small and noisy. In the estimates from Gs, the localized response is recovered however the remaining component is contains waves because the BOLD response does where as Gs does not. Hence, this leads to a larger and smooth time to peak profile that has wave properties.

Summary and discussion Summary In this paper, we have developed a method to infer neuronal activity through deconvolution of surface-reconstructed fMRI data using spatiotemporal hemodynamic response functions (stHRFs). This paper focuses on using an stHRF obtained from a first principles physiological model of hemodynamics — Gp, and compared theses inferences to those obtained using with the traditional space time separable stHRF Gs. We first used simulated BOLD data generated by Gp to validate and analyze our spatiotemporal deconvolution method on simulated data that contained damped traveling waves. Deconvolving these test BOLD data using Gp yielded estimates that closely approximated the input drive. Upon deconvolution with the proposed space–time separable stHRF, Gs i.e., the mismatched filter, yielded false neural dynamics that included spurious slow neural waves and “ghost” neural drives. This finding stresses the importance of using the correct stHRF, i.e., one based on the underlying physiology. We then applied these deconvolution methods to empirical BOLD data acquired from the primary visual cortex during a stimulusevoked paradigm. Previous analyses of these data in surface-based reconstructions revealed slow traveling hemodynamic waves (Aquino

K.M. Aquino et al. / NeuroImage 94 (2014) 203–215

(a)

ζ ⊗ GP

(b)

ζ

211

20

20

0.8

0.8 15

15

0.6

10

0.4

5

t (s)

t (s)

0.6

0.2

0 −5

0 x (mm)

0.4

5

0.2 0

0 −5

5

ζ using GP

(c)

10

0 x (mm)

ζ using GS

(d)

20

20

0.8

0.8 0.6

0.6

15

5

15

0.4

10

0.2

5

−0.2 0 x (mm)

0.2

10

0 −0.2

5

0

0 −5

t (s)

t (s)

0.4

5

−0.4 0 −5

0 x (mm)

5

Fig. 11. The deconvolution of the BOLD response to two neural sources. (a) The neural drive using Eq. (19). (b) The predicted bold response by convolving ζ with Gp. (c) The estimate of ζ^ using Gp. (d) The estimate of ζ^ using Gs.

et al., 2012). Deconvolution with Gp yielded a local, spatially stationary neural drive, consistent with independent recordings. Two points are noteworthy here: Firstly, we deconvolve the wave components of the BOLD signal to infer a nonpropagating neural source. Secondly, the spatial extent of the inferred neural drive is narrower than that of the BOLD response. Hence the deconvolution is similar in effect to the “sharpening effect” now applied almost universally to diffusion imaging data (Tournier et al., 2007). Notably, our deconvolution method is derived from the physiology. Crucially, where BOLD waves are present in highly resolved data (such as in our experiment) deconvolution with the unphysical Gs would imply existence of slow neural source waves traveling at speeds of the order of 1 mm s−1, which are not supported by independent recordings. If waves are not present in the empirical data (or are weak, such as in some of our subjects) our approach would still function to more accurately estimate the underlying neuronal response [see Fig. 2 of Aquino et al. (2012)]. These findings are complementary to a recent study by Bießmann et al. (2012) which showed higher accuracy when using a nonseparable stHRF (albeit not based directly on the physiology) compared to a separable one. Hence, we conclude that using Gp will estimate neural activity more accurately.

0.4

20

15

0.2

15

10

0

5 0 −5

0 x (mm)

5

In the present study we modeled and deconvolved the BOLD response using a linear approximation. Here, we address the validity of this approximation. It has been previously demonstrated theoretically that that as long as the neural drive is sufficiently small or modest, the BOLD response is approximately linear (Aquino et al., 2012; Drysdale et al., 2010; Robinson et al., 2006). Therefore the linear convolution model to determine the BOLD response is restricted to low amplitude responses, i.e., neuronal responses likely to hold in ecologically likely scenarios (Attwell and Iadecola, 2002; Schölvinck et al., 2008). Hence this method does not take into account effects such as a nonlinear relationship between stimulus strength and response amplitude or delayed compliance of cortical vessels for long stimulations (Buxton et al., 2004; Mandeville et al., 1999). The effects that these will have on the neural estimates require future work. However, for the current evoked responses, theoretical analyses suggest that these nonlinearities are unlikely to cause large differences in the estimated effects (Aquino et al. 2014).

(c)

ζ using GP 0.5

0

10

−0.2

5

−0.4

0 −5

ζ using GS

20

−0.5 0 x (mm)

5

0.5

15 t (s)

(b)

BOLD

20

t (s)

t (s)

(a)

Linear approximations

0

10 5 0 −5

−0.5 0 x (mm)

5

Fig. 12. The deconvolution of evoked fMRI data using Gp and Gs. (a) The evoked spatiotemporal BOLD response of a single subject projected onto one spatial dimension. (b) The deconvolved estimate using Gp. (c) The deconvolved estimate using Gs. Each response has been normalized to have a maximum amplitude of unity, and hence dimensionless.

212

K.M. Aquino et al. / NeuroImage 94 (2014) 203–215

(b) 1

10

0.8

8

time to peak (s)

max amplitude

(a)

0.6 0.4 0.2 0 −5

BOLD

6

GP

4

GS

2

0

0 −5

5

x (mm)

0

5

x (mm)

Fig. 13. Similarities and differences between the deconvolved estimates ζ^ . (a) Shows the maximum amplitude as a function of distance, and (b) shows the time to peak as a function of distance. The estimates using Gp and Gs are in black and blue respectively, and the red point set is the BOLD response.

Implications for existing studies

Gp vs Gs

The effects presented in this study become pronounced and significant at highly resolved fMRI data. Firstly, the presented hemodynamic waves are damped and therefore have a limited spatial scale. From independent estimates of the model parameters (Aquino et al., 2012, 2014) we show that the waves can travel up to 5–10 mm before they decay and can induce signal delays up to 5 s. Therefore, the difference between deconvolution estimates from the physiological to the assumed stHRF will thus affect inversions of neural activity at spatial scales in the order of millimeters. Therefore, it is unlikely that inferences of false neural waves and ghost neural responses will be present in current causal modeling techniques since they compare regions that are at least centimeters apart. However, as fMRI resolutions advance, data will contain rich spatiotemporal dynamics at the scale of b1 mm and b 1 s, and from our work we show that as it stands, direct application of current casual modeling to such high resolution data would produce false inferences, essentially counteracting improvements of resolution. Although our study develops a method to deconvolve delays due to hemodynamic waves, we stress that this technique does not remove all sources of hemodynamic delay. Some of the delays not covered include the anti correlated BOLD signal due to large draining veins (Lee et al., 1995), the delays between hemodynamic layers (Tian et al., 2010), and any apparent delay caused by the aliasing of high frequency components onto the BOLD response (Kriegeskorte et al., 2010). These delays will have to be deconvolved at a different stage and cause distinctly different properties on the BOLD response.

Our study shows that the dynamics of the neural estimate depends on the choice of the HRF, here we detail the strengths and weakness of Gp and Gs. Firstly, in regards to parameters: Gp is derived from a physical model containing 13 independent parameters, all of which have a physiological origin and are therefore not free. In contrast Gs contains 10 parameters one of which, the spatial width of the point spread function, is a free parameter. We additionally showed in our simulations that Gp still outperforms Gs with this increased number of parameters, however independent recordings of the neural drive are needed to test this in empirical data. Secondly, in regards to artifacts, we showed in the simulations that if the mismatched filter Gs is used to deconvolve data with waves, this will imply the existence of false neural waves and large ghost responses. For completeness, if we assume that Gp is the mismatched filter, the neural estimates will appear smoothed in time and in space, as described in Appendix B. That is the precision of the inferences will be diminished (in space) but there will not be a systematic bias. Thirdly, an important point is regarding the differences of the time to peak profiles between estimates using Gp vs Gs. Although these are small in our simulations and data, their differences imply distinctly different mechanisms. In particular, deconvolution with Gs on real data reveals a time to peak profile that would imply neural waves traveling with velocity ≈ 1 mm s−1, which is unsupported by independent recordings. This is at least 1–2 orders of magnitude slower than neural waves measured by independent approaches (Benucci et al., 2007; Grinvald et al., 1994; Xu et al., 2007), in contrast to the estimates using Gp, which imply a localized response consistent with independent recordings.

(a)

(b) 8

1

time to peak (s)

max amplitude

0.8 0.6 0.4 0.2 0 −5

0

x (mm)

5

6

BOLD 4

GP GS

2

0 −5

0

5

x (mm)

Fig. 14. The similarities and differences of the deconvolution estimates calculated for the entire data set. (a) Is the maximum amplitude as a function of distance and (b) is the time to peak as a function of distance averaged across all subjects. The estimates using Gp and Gs are in black and blue respectively, the red point set is the BOLD response and the error bars indicate S.E.M. The fits for the time to peak, TTP as a function of x were given by TTP (x) = −2.4, and TTP (x) = 0.088x − 0.010 + 2.4(R2 = 0.91) for the estimates from Gp and Gs respectively.

K.M. Aquino et al. / NeuroImage 94 (2014) 203–215

213

and active neural contributions to flow control [(Devor et al., 2008), see Attwell et al. (2010) for a review]. Such innovations would enable one to determine neural firing as opposed to the overall neural drive. In addition this could be addressed by using a neural field model combined with EEG–fMRI to refine estimates of neural activity (Riera et al., 2006; Rosa et al., 2011). Although there will ultimately be a need to go beyond the present work by taking into account individual local vascular geometry, the present work provides a method to probe spatiotemporal neural activity at available fMRI resolutions.

Finally, the use of Gs is further supported by prior analysis of high resolution fMRI data. In particular, Shmuel et al. (2007) showed that although a separable approximation (i.e., Gs) could account for a high proportion of the total variance in a measured stHRF, the spatial width of the stHRF increased with time as would be expected when propagating waves are present: This supports our conclusion that the stHRF is not separable. Indeed, it is only at late post-stimulus times that one would expect non-separability to be clearly observable. The total variance is not sensitive to this aspect. To conclude, we note that although we have separately defined Gs, this HRF can be derived directly from the spatiotemporal hemodynamic model as a limiting case by having a combination of having a high hemodynamic velocity vβ and high damping coefficient Γ as seen in (Aquino et al., 2012, 2014). For example, the post-stimulus undershoot seen in some empirical data and with our simulated Gs can also be simulated with Gp, using appropriate parameters. Hence all effects are covered by using Gp.

Conclusions As advances in fMRI scanner hardware continue to improve the spatiotemporal resolution of the BOLD signal, systematic analyses of the spatiotemporal properties will become critical to interpret BOLD activity and extract maximum information and resolution from fMRI. In the first instance, our approach would likely find most relevance to local surface-reconstructed data acquired with high spatial resolution protocols. For example, several groups have used fMRI to unravel finescale spatiotemporal patterns of neural activity; e.g., in the detection of orientation in primary visual cortex (Haynes and Rees, 2005; Kamitani and Tong, 2005; Swisher et al., 2010). As a “sharpening deconvolution”, our technique would be able to refine such maps, separating neural and hemodynamic contributions to these measures. In the longer term, an accurate model of spatiotemporal hemodynamics has far broader application. For example, a forward model is crucial to data inversion schemes such as the approach implemented in DCM (Friston et al., 2003). In concert with neural field models of brain activity (Deco et al., 2008), the present stHRF will enable new model-driven inferences of neuronal activity within cortical regions, and of interactions between regions.

Improvements to the deconvolution There are a number of possible improvements to the proposed deconvolution technique. First, a more systematic approach can be used to describe the NSR spectrum either by a quantitative model [e.g. Kerr et al. (2009)] or with tighter constraints on the G spectrum using independent observations from the same subject. Second, as described by Gitelman et al. (2003), the parameter estimates could be optimized so that the difference between the predicted BOLD signal Y, and the predicted neural estimate ζ^ with the stHRF G is minimized    i.e., by minimizing Y−ζ^ ⊗G. This requires a model for the spatiotemporal neural response, and the results presented in this study will aid in determining a model for this response. Refinements to the deconvolution can also be made via including additional effects on modeling, in addition to the discussion on nonlinearities above another point is to increase the accuracy of the neurovascular coupling. These additions include specifically taking into account astrocytic delays (Araque et al., 2001; Carmignoto and Gómez-Gonzalo, 2010; Haydon and Carmignoto, 2006; Zonta et al., 2003) as well as direct

Acknowledgments The authors thank Cliff Kerr for fruitful discussions on deconvolution. The Australian Research Council and the Westmead Millennium Institute supported this work.

Appendix A. Table of model parameters This appendix details the model parameters of the spatiotemporal hemodynamic response functions and the nominal parameters used for the simulations.

Table A.1 The model variables and parameters. In each row, the five columns detail the quantity, its symbolic formula, its nominal range from physiology, the units, and the source of the parameter range, respectively. The last column details notes on each variable/parameter: (a) denotes a quantity that does not have a direct analog in the Balloon Model, (b) denotes a parameter used in previous Balloon Models, and (c) denotes the set of independent variables needed to calculate the stHRF. Parameter

Symbol

Value/range

Units

Source

Notes

Blood flow signal decay rate Flow-dependent elimination constant Grubb exponent Mean elasticity exponent of cortical vessels Hemodynamic transit time Resting blood oxygen extraction fraction Oxygen consumption rate Ratio of hemoglobin concentration to blood density Resting blood volume fraction Magnetic field parameters in the signal equation at 3 T, TE = 30 ms Natural frequency of flow response Wave propagation velocity Effective blood viscosity Wave damping rate Cortical thickness Perpendicular spatial frequency

0.65 0.41 0.31 3.2 1–4 0.4 0.4 0.15 0.03 4.2,1.7,0.41 0.56 1–20 106–850 0.2 − 2 ≈3 ≈200

s−1 s−2 – – s s s−1 mol kg−1 – – s−1 mm s−1 kg m−3 s−1 s−1 mm m−1

Friston et al. (2003) Friston et al. (2003) Friston et al. (2003) Drysdale et al. (2010) Buxton et al. (2004) Buxton et al. (2004) Aquino et al. (2012) Aquino et al. (2012) Friston et al. (2003) Huppert et al. (2006) Aquino et al. (2012) Aquino et al. (2012) Aquino et al., (2014) Aquino et al. (2012) Duvernoy et al. (1981) Aquino et al., (2014)

(b) (b) (b) (a) (b) (b) (b),(c) (a) (b),(c) (b),(c) (b),(c) (a),(c) (a) (a),(c) (a) (a),(c)

Effective spatial frequency

κ γ α β = 1/α τ ρHb η = ρHb/τ ψ V0 k1,k2,k3 ωf qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi β−1 ν β ¼ c1 c2 ξ β D
 Γ ¼ 12 ρD þ βτ L k0 = cos−1(0.8)/L qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 kz ¼ k0 þ ν1 C Z βτρD

≈200–950

m−1

Aquino et al., (2014)

(a),(c)

Outflow normalization constant

C Z ¼ 1=k1mm sinðk LÞ

≈0.36

–

Aquino et al., (2014)

(a),(c)

f

2 β

0

0

f

214

K.M. Aquino et al. / NeuroImage 94 (2014) 203–215

(b) 0.5 0 −0.5

10

−1

5

(d)

20 0.5

15

0.6

15

0.8

20

0.8

0.6

15

0.4

t (s)

t (s)

15

(c)

20

10 0 5

−1.5

t (s)

20

10

0.2 0

5

t (s)

(a)

0.4 10 0.2 5

0

−0.2 −2

0 −5

0

0 −5

5

0

x (mm)

−0.5

5

x (mm)

0

−5

0

x (mm)

5

0 −5

0

5

−0.2

x (mm)

Fig. 15. Estimates of the neural drive by creating a HRF using Gs and then deconvolving with a mismatched filter Gp. Panels (a)–(d) are deconvolved estimates of ζ^ using (a) fc = 2 Hz, kc/(2π) = 1000 m−1 (b) fc = 0.1 Hz, kc/(2π) = 500 m−1, (c) fc = 0.1 Hz, kc/(2π) = 200 m−1, and (d) fc = 0.05 Hz, kc/(2π) = 200 m−1. Each response has been normalized to have a maximum amplitude of unity, and hence dimensionless.

(a)

(b)

20

(c)

20

20

0.8

0.8 15

10

0.4 0.2

5

10

0.4 0.2

5

0 0 −10

0

x (mm)

10

0.8 15

0.6

t (s)

0.6

t (s)

t (s)

15

0.6

10

0.4 0.2

5

0 0 −10

0

x (mm)

10

0 0 −10

0

10

x (mm)

Fig. 16. The effect of neural width on synthetic BOLD responses. The BOLD responses are simulated using Eq. (18) using widths that have FWHM of (a) 1 mm, (b) 4 mm, (c) 6 mm. Each of these responses are normalized to have maximum amplitude of unity, hence dimensionless.

Appendix B. Deconvolution using Gp as the mismatched filter The purpose of this appendix is to show the artifacts that arises if Gp is the mismatched filter in the deconvolution. For completeness, we repeated the analysis in Deconvolution of synthetic data by using synthetic data generated by Gs and deconvolving the data with Gs. Fig. 15 shows the result of this analysis. It is apparent visually, that the neural drive is not accurately estimated [c.f. Fig. 3(b)] and contains artifacts in its reconstruction. The key artifact in these estimates is the overestimation of the spatial width of ζ. The reason why this occurs is because Gp contains both local and propagating modes (Aquino et al. 2014). By increasing the spatial extent of the neural drive, the local nonpropagating mode can dominate the response, and as shown in Fig. 16 the responses have less propagating components for this simple neural drive. Hence deconvoluting a non-propagating BOLD response with Gp, the spatial extent of ζ^ will be over estimated. Appendix C. Supplementary data Supplementary data to this article can be found online at http://dx. doi.org/10.1016/j.neuroimage.2014.03.001. References Akaike, H., 1974. A new look at the statistical model identification. IEEE Trans. Autom. Control 19, 716–723. Aquino, K.M., Schira, M.M., Robinson, P.A., Drysdale, P.M., Breakspear, M., 2012. Hemodynamic traveling waves in human visual cortex. PLoS Comput. Biol. 8, e1002435. Aquino, K.M., Robinson, P.A., Drysdale, P.M., 2014. Spatiotemporal hemodynamic response functions derived from physiology. J. Theor. Biol. 347, 118–136. Araque, A., Carmignoto, G., Haydon, P.G., 2001. Dynamic signaling between astrocytes and neurons. Annu. Rev. Physiol. 63, 795–813. Attwell, D., Iadecola, C., 2002. The neural basis of functional brain imaging signals. Trends Neurosci. 25, 621–625. Attwell, D., Buchan, A.M., Charpak, S., Lauritzen, M., MacVicar, B.A., Newman, E.A., 2010. Glial and neuronal control of brain blood flow. Nature 468, 232–241.

Benucci, A., Frazor, R.A., Carandini, M., 2007. Standing waves and traveling waves distinguish two circuits in visual cortex. Neuron 55, 103–117. Bianciardi, M., Fukunaga, M., Van Gelderen, P., de Zwart, J.A., Duyn, J.H., 2010. Negative BOLD-fMRI signals in large cerebral veins. J. Cereb. Blood Flow Metab. 31, 401–412. Bießmann, F., Murayama, Y., Logothetis, N., Müller, K., Meinecke, F., 2012. Improved decoding of neural activity from fMRI signals using non-separable spatiotemporal deconvolutions. Neuroimage 61, 1031–1042. Boynton, G.M., Engel, S.A., Glover, G.H., Heeger, D.J., 1996. Linear systems analysis of functional magnetic resonance imaging in human V1. J. Neurosci. 16, 4207–4221. Buxton, R.B., Wong, E.C., Frank, L.R., 1998. Dynamics of blood flow and oxygenation changes during brain activation: the balloon model. Magn. Reson. Med. 39, 855–864. Buxton, R.B., Uludaug, K., Dubowitz, D.J., Liu, T.T., 2004. Modeling the hemodynamic response to brain activation. Neuroimage 23, S220–S233. Carmignoto, G., Gómez-Gonzalo, M., 2010. The contribution of astrocyte signalling to neurovascular coupling. Brain Res. Rev. 63, 138–148. Chang, C., Thomason, M.E., Glover, G.H., 2008. Mapping and correction of vascular hemodynamic latency in the BOLD signal. Neuroimage 43, 90–102. Chen, Y., Namburi, P., Elliott, L., Heinzle, J., Soon, C., Chee, M., Haynes, J., 2011. Cortical surface-based searchlight decoding. Neuroimage 56, 582–592. David, O., Guillemain, I., Saillet, S., Reyt, S., Deransart, C., et al., 2008. Identifying neural drivers with functional MRI: an electrophysiological validation. PLoS Biol. 6, 2683–2697. Deco, G., Jirsa, V., Robinson, P., Breakspear, M., Friston, K., 2008. The dynamic brain: from spiking neurons to neural masses and cortical fields. PLoS Comput. Biol. 4, e1000092. Devor, A., Hillman, E.M., Tian, P., Waeber, C., Teng, I.C., Ruvinskaya, L., Shalinsky, M.H., Zhu, H., Haslinger, R.H., Narayanan, S.N., et al., 2008. Stimulus-induced changes in blood flow and 2-deoxyglucose uptake dissociate in ipsilateral somatosensory cortex. J. Neurosci. 28, 14347–14357. Drysdale, P.M., Huber, J., Robinson, P.A., Aquino, K.M., 2010. Spatiotemporal BOLD dynamics from a porous elastic hemodynamic model. J. Theor. Biol. 265, 524–534. Duvernoy, H.M., Delon, S., Vannson, J.L., 1981. Cortical blood vessels of the human brain. Brain Res. Bull. 7, 519–579. Engel, S.A., Glover, G.H., Wandell, B.A., 1997. Retinotopic organization in human visual cortex and the spatial precision of functional MRI. Cereb. Cortex 7, 181–192. Frank, L.R., Buxton, R.B., Wong, E.C., 2000. What's in the noise now? Int. Soc. Magn. Reson. Med. 8–812. Frank, L.R., Buxton, R.B., Wong, E.C., 2001. Estimation of respiration-induced noise fluctuations from undersampled multislice fMRI data. Magn. Reson. Med. 45, 635–644. Friston, K.J., 1995. Regulation of rCBF by diffusible signals: an analysis of constraints on diffusion and elimination. Hum. Brain Mapp. 3, 56–65. Friston, K.J., 2009. Causal modelling and brain connectivity in functional magnetic resonance imaging. PLoS Biol. 7, e1000033. 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, 466–477. Friston, K.J., Harrison, L., Penny, W., 2003. Dynamic causal modelling. Neuroimage 19, 1273–1302. Gardner, J.L., 2010. Is cortical vasculature functionally organized? Neuroimage 49, 1953–1956. Gitelman, D.R., Penny, W.D., Ashburner, J., Friston, K.J., 2003. Modeling regional and pyschophysiologic interactions in fMRI: the importance of hemodynamic deconvolution. Neuroimage 19, 200–207. Glover, G.H., 1999. Deconvolution of impulse response in event-related BOLD fMRI. Neuroimage 9, 416–429. Grinvald, A., Lieke, E.E., Frostig, R.D., Hildesheim, R., 1994. Cortical point-spread function and long-range lateral interactions revealed by real-time optical imaging of macaque monkey primary visual cortex activity. J. Neurosci. 14, 2545–2568. Handwerker, D.A., Ollinger, J.M., D'Esposito, M., 2004. Variation of bold hemodynamic responses across subjects and brain regions and their effects on statistical analyses. Neuroimage 21, 1639–1651. Haydon, P.G., Carmignoto, G., 2006. Astrocyte control of synaptic transmission and neurovascular coupling. Physiol. Rev. 1009–1031. Haynes, J.D., Rees, G., 2005. Predicting the orientation of invisible stimuli from activity in human primary visual cortex. Nat. Neurosci. 8, 686–691. Horwitz, B., 2003. The elusive concept of brain connectivity. Neuroimage 19, 466–470.

K.M. Aquino et al. / NeuroImage 94 (2014) 203–215 Huppert, T.J., Hoge, R.D., Diamond, S.G., Franceschini, M.A., Boas, D.A., 2006. A temporal comparison of BOLD, ASL, and NIRS hemodynamic responses to motor stimuli in adult humans. Neuroimage 29, 368–382. Kamitani, Y., Tong, F., 2005. Decoding the visual and subjective contents of the human brain. Nat. Neurosci. 8, 679–685. Kerr, C.C., Rennie, C.J., Robinson, P.A., 2009. Deconvolution analysis of target evoked potentials. J. Neurosci. Methods 179, 101–110. Kriegeskorte, N., Cusack, R., Bandettini, P., 2010. How does an fMRI voxel sample the neuronal activity pattern: compact-kernel or complex spatiotemporal filter? Neuroimage 49, 1965–1976. Krueger, G., Glover, G.H., 2001. Physiological noise in oxygenation-sensitive magnetic resonance imaging. Magn. Reson. Med. 46, 631–637. Lee, A.T., Glover, G.H., Meyer, C.H., 1995. Discrimination of large venous vessels in timecourse spiral blood-oxygen-level-dependent magnetic-resonance functional neuroimaging. Magn. Reson. Med. 33, 745–754. Mandeville, J.B., Marota, J.J.A., Ayata, C., Zaharchuk, G., Moskowitz, M.A., Rosen, B.R., Weisskoff, R.M., 1999. Evidence of a cerebrovascular postarteriole windkessel with delayed compliance. J. Cereb. Blood Flow Metab. 19, 679–689. Melrose, D.B., McPhedran, R.C., 1991. Electromagnetic Processes in Dispersive Media: A Treatment Based on the Dielectric Tensor. Cambridge Univ. Press, Cambridge. Menon, R.S., Goodyear, B.G., 1999. Submillimeter functional localization in human striate cortex using BOLD contrast at 4 Tesla: implications for the vascular point-spread function. Magn. Reson. Med. 41, 230–235. Parkes, L.M., Schwarzbach, J.V., Bouts, A.A., Deckers, R.h.R., Pullens, P., Kerskens, C.M., Norris, D.G., 2005. Quantifying the spatial resolution of the gradient echo and spin echo BOLD response at 3 Tesla. Magn. Reson. Med. 54, 1465–1472. Riera, J.J., Wan, X., Jimenez, J.C., Kawashima, R., 2006. Nonlinear local electrovascular coupling. I: A theoretical model. Hum. Brain Mapp. 27, 896–914. Robinson, P.A., Drysdale, P.M., Van der Merwe, H., Kyriakou, E., Rigozzi, M.K., Germanoska, B., Rennie, C.J., 2006. BOLD responses to stimuli: dependence on frequency, stimulus form, amplitude, and repetition rate. Neuroimage 31, 585–599. Rosa, M.J., Kilner, J.M., Penny, W.D., 2011. Bayesian comparison of neurovascular coupling models using EEG–fMRI. PLoS Comput. Biol. 7, e1002070.

215

Schölvinck, M.L., Howarth, C., Attwell, D., 2008. The cortical energy needed for conscious perception. Neuroimage 40, 1460–1468. Schwarz, G., 1978. Estimating the dimension of a model. Ann. Stat. 6, 461–464. Shmuel, A., Yacoub, E., Chaimow, D., Logothetis, N.K., Ugurbil, K., 2007. Spatio-temporal pointspread function of fMRI signal in human gray matter at 7 Tesla. Neuroimage 35, 539–552. Stephan, K.E., Roebroeck, A., 2012. A short history of causal modelling of fMRI data. Neuroimage 62, 856–863. Swisher, J.D., Gatenby, J.C., Gore, J.C., Wolfe, B.A., Moon, C.H., Kim, S.G., Tong, F., 2010. Multiscale pattern analysis of orientation-selective activity in the primary visual cortex. J. Neurosci. 30, 325–330. Tian, P., Teng, I.C., May, L.D., Kurz, R., Lu, K., Scadeng, M., Hillman, E., De Crespigny, A.J., D'Arceuil, H.E., Mandeville, J.B., Marota, J.J.A., Rosen, B., Liu, T.T., Boas, D.A., Buxton, R.B., Dale, A.M., Devor, A., 2010. Cortical depth-specific microvascular dilation underlies laminar differences in blood oxygenation level-dependent functional MRI signal. Proc. Natl. Acad. Sci. U. S. A. 107, 15246–15251. Tournier, J., Calamante, F., Connelly, A., et al., 2007. Robust determination of the fibre orientation distribution in diffusion MRI: non-negativity constrained super-resolved spherical deconvolution. Neuroimage 35, 1459–1472. Vanduffel, W., Tootell, R.B., Schoups, A.A., Orban, G.A., 2002. The organization of orientation selectivity throughout macaque visual cortex. Cereb. Cortex 12, 647–662. Wang, S.J., Luo, L.M., Liang, X.Y., Gui, Z.G., Chen, C.X., 2005. Estimation and removal of physiological noise from undersampled multi-slice fMRI data in image space. Conf. Proc. IEEE Eng. Med. Biol. Soc. 1371–1373. Worsley, K.J., Friston, K.J., 1995. Analysis of fMRI time-series revisited—again. Neuroimage 2, 173–181. Xu, W., Huang, X., Takagaki, K., Wu, J.Y., 2007. Compression and reflection of visually evoked cortical waves. Neuron 55, 119–129. Zarahn, E., Aguirre, G.K., D'Esposito, M., 1997. Empirical analyses of BOLD fMRI statistics. Neuroimage 5, 179–197. Zonta, M., Angulo, M., Gobbo, S., Rosengarten, B., Hossmann, K.A., Pozzan, T., Carmignoto, G., 2003. Neuron-to-astrocyte signalling is central to the dynamic control of brain microcirculation. Nat. Neurosci. 6, 43–50.