Home | About | Journals | Submit | Contact Us | Français |

**|**J R Soc Interface**|**v.13(118); 2016 May**|**PMC4892270

Formats

Article sections

- Abstract
- 1. Introduction
- 2. Spatio-temporal haemodynamic model
- 3. Natural modes of the blood oxygen-level dependent response
- 4. Features of components Tj(k,ω)
- 5. Blood oxygen-level dependent response to a localized neural activity: Green function
- 6. Illustrative applications
- 7. Summary and conclusions
- References

Authors

Related links

J R Soc Interface. 2016 May; 13(118): 20160253.

PMCID: PMC4892270

e-mail: ua.ude.yendys@gnap.semaj

Received 2016 March 31; Accepted 2016 April 12.

Copyright © 2016 The Author(s)

Published by the Royal Society. All rights reserved.

The blood oxygen-level dependent (BOLD) response to a neural stimulus is analysed using the transfer function derived from a physiologically based poroelastic model of cortical tissue. The transfer function is decomposed into components that correspond to distinct poles, each related to a response mode with a natural frequency and dispersion relation; together these yield the total BOLD response. The properties of the decomposed components provide a deeper understanding of the nature of the BOLD response, via the components' frequency dependences, spatial and temporal power spectra, and resonances. The transfer function components are then used to separate the BOLD response to a localized impulse stimulus, termed the Green function or spatio-temporal haemodynamic response function, into component responses that are explicitly related to underlying physiological quantities. The analytical results also provide a quantitative tool to calculate the linear BOLD response to an arbitrary neural drive, which is faster to implement than direct Fourier transform methods. The results of this study can be used to interpret functional magnetic resonance imaging data in new ways based on physiology, to enhance deconvolution methods and to design experimental protocols that can selectively enhance or suppress particular responses, to probe specific physiological phenomena.

Functional magnetic resonance imaging (fMRI), based on the blood oxygen-level dependent (BOLD) signal, has become increasingly popular in the last decades because of its ability to probe the function of various regions of the brain. The experiments are performed through non-invasive measurement of vascular oxygenation changes in response to neural activity [1–5]. The dynamics of the BOLD response can be further understood by formulating a haemodynamic response function (HRF) that represents the spatio-temporal properties of the response to a localized impulse of stimulation; for neuroimaging, the stimulation is related to the neural activity. Such a formulation will potentially enable better mapping of the brain areas that respond to specific experimental stimuli.

Attempts have been made to model the BOLD signal by linking changes in haemodynamic quantities, such as cerebral blood flow (CBF), cerebral blood volume (CBV) and deoxygenated haemoglobin (dHb). One of them is the temporal ‘balloon model’ [6] that treats the cerebral vasculature as an elastic balloon that inflates/deflates during an increase/decrease in blood volume. Paired with the neurovascular model of [7], this has successfully predicted observed temporal properties of the BOLD signal, such as stimulus-dependent increases and undershoots [7,8]. Many studies on the spatio-temporal variation of BOLD, on the other hand, have rested on an assumption that the spatio-temporal haemodynamic response function (stHRF) is just a product of a temporal HRF and a spatial point spread function. The spatial function has been widely approximated by a Gaussian function [9], motivated by numerical fitting in three-dimensional fMRI experiments [10,11] and rapid spatial diffusion of vasodilatory agents, such as nitric oxide [12]. However, the simple product of temporal and spatial HRFs ignores intrinsic brain dynamics such as blood flow in between neighbouring veins and arteries, finite blood propagation speeds and haemodynamic waves [13].

A recent physiologically based model, which treats cortical tissue as a poroelastic medium, is successful in predicting both the spatial and temporal characteristics of the BOLD response [13–15]. The model was formulated using the methods from the theory of linear poroelasticity in geophysics [16], resulting in a transfer function that predicts the linear BOLD response to a neural activity. It comprises coupled nonlinear partial differential equations that embody the mean field approximations of haemodynamic changes of blood flow, blood mass, pore pressure and haemoglobin concentration in response to neural activity as governed by physical conservation laws.

The features of the transfer function derived from poroelastic haemodynamics relate the dynamics of inputs to determine the spatio-temporal properties of the BOLD response. One important feature is the frequency poles that allow recovery of the natural frequencies that govern the inherent modes of oscillation of the BOLD signal. In order to find the poles, one can use the partial fraction method widely used in control systems theory [17] and computer modelling [18]. The technique transforms a transfer function into linear combinations of lower-order components, thereby reducing the complexity of the system under study. The objective of this work is to analytically decompose the BOLD transfer function into partial fractions to derive the poles and response modes. This may enable us to formulate optimization protocols that allow one to selectively enhance or suppress specific modes/responses depending on which physiological phenomenon is relevant.

The paper is organized as follows: in §2, we summarize the spatio-temporal haemodynamic model [13–15] and discuss how the resulting transfer function predicts the linear BOLD response to a neural activity. In §3, we detail the partial fraction method to decompose the BOLD response into natural modes, each corresponding to a frequency pole. In §4, we characterize the spatio-temporal frequency properties of these components in terms of frequency responses, power spectra and resonances. In §5, we derive analytic expressions for the BOLD response to localized impulsive neural activity, mathematically modelled as a Dirac delta function; we term this response the Green function or the stHRF. The results provide a powerful quantitative tool to calculate the linear BOLD response to an arbitrary neural drive, which is faster and more efficient than traditional Fourier transform methods. In §6, we apply our analytical results to determine the BOLD response to a spatio-temporal Gaussian stimulus and compare with experimental fMRI data. We show how the response can be decomposed into its components. We also show an illustrative example to selectively enhance the BOLD response components by using a spatio-temporally moving sinusoidal neural stimulus. Finally, in §7, we summarize our major results and provide key points in designing appropriate protocols to test our findings in future experiments.

Many studies have shown that the BOLD signal, measured using fMRI, can serve as a proxy for neural responses [4,5]. The physical principles underlying fMRI experimental data can be better understood by formulating a theory that can predict the haemodynamic response to neural activity. In this section, we outline the key results of a recent haemodynamic model derived from physiological properties of cortical tissue that produces a transfer function that can predict the linear BOLD response to a neural activity. Readers are referred to the original references for detailed description and derivation [13–15].

Here we outline and briefly discuss the spatio-temporal haemodynamic model. Its equations can be used to predict the spatio-temporal BOLD signal *Y*.

The model starts by considering a cortical tissue approximated as a poroelastic medium in which the vessels serve as the ‘pores’ of the system as shown in figure 1. In the presence of a neural activity *ζ*(**R**, *t*) at location **R** at time *t*, the neurovascular coupling between neurons and blood vessels leads to changes in blood inflow. The arterial blood inflow rate *F*(**R**, *t*) is modelled as a damped harmonic oscillator driven by the neural drive *ζ*(**R**, *t*) [7,14]

2.1

where *κ*, *γ* and *F*_{0} are the blood flow signal decay rate, the flow-dependent elimination constant and the resting blood inflow, respectively.

Arbitrary volume element of poroelastic cortical tissue. Blood flows through the vessels, serving as pores of the tissue, in the direction of the solid arrows. Broken arrow indicates the average direction of the flow with velocity **v**. (Online version in **...**

The response owing to inflows or outflows is constrained by changes in blood vessel pressure and conservation of mass and momentum of blood within the tissue. This is collectively modelled by the nonlinear haemodynamic wave equation for blood mass Ξ(**R**, *t*) [13,15]

2.2

where *D* quantifies damping owing to effective blood viscosity, *ρ _{f}* is the blood density,

2.3

where *c*_{2} is the porous coupling constant and *β* is the elasticity exponent of cortical vessels. In this paper, we associate *β* with the reciprocal of Grubb's exponent, which relates CBF to CBV, such that *β* ≈ 3.23 [14,19]. This *β* value implies that the average behaviour of the vessels is hyperelastic (i.e. more resistant to stretching than perfectly elastic vessels).

Analogous to the conservation of blood mass in the arterioles, conservation of dHb must also be followed in the veins. Changes in the concentration of dHb *Q*(**R**, *t*) are described by the dynamical equation

2.4

where **v**(**R**, *t*) is the average blood velocity, **v*** _{F}*(

2.5

and

2.6

whereas conservation of momentum governs the average blood velocity such that

2.7

Finally, the BOLD signal *Y* is modelled by the semi-empirical relation [20]

2.8

where *V*_{0} is the resting blood volume fraction, *Q*_{0} is the resting dHb concentration per unit cerebral cortical volume, and *k*_{1}–*k*_{3} are parameters dependent on the configuration of the fMRI scanner used. Equation (2.8) expresses the spatio-temporal BOLD signal *Y*(**R**, *t*) in terms of the quantities *F*(**R**, *t*), Ξ(**R**, *t*), *P*(**R**, *t*) and *Q*(**R**, *t*) for a given neural drive *ζ*(**R**, *t*).

Equations (2.1)–(2.8) can be used to calculate the nonlinear haemodynamic response to any kind of neural drive *ζ*(**R**, *t*). However, when the neural drive is short in duration and low in amplitude, it has been shown in experiments that the haemodynamic response is approximately linear [13,21]. In this case, we can treat the haemodynamic variables *F*, Ξ, *Q*, **v**, *P*, *ζ* and *Y* as being linear perturbations from their steady-state values allowing us to linearize equations (2.1)–(2.8).

Fourier methods can be used to understand the properties of the linearized equations at a given spatial frequency **k** and temporal angular frequency *ω*. This leads to a set of transfer functions *T _{θζ}* which describes the change of the haemodynamic quantity

2.9

where *Y*(**k**, *ω*) and *ζ*(**k**, *ω*) are the Fourier transforms of the BOLD signal *Y*(**R**, *t*) and neural drive *ζ*(**R**, *t*), respectively.

Upon averaging of dynamics through the thickness of the cortex, linearization of equations (2.1)–(2.8), and mapping to Fourier space, equation (2.9) becomes [15]

2.10

where is the spatial wavenumber in the plane of the cortex, *C _{z}* is the outflow normalization constant,

Equations (2.9) and (2.10) can be used to calculate the BOLD response to an arbitrary drive. This is achieved by taking the inverse Fourier transform of *Y*(**k**, *ω*) to obtain the two-dimensional linear BOLD response in coordinate space as

2.11

where from here on, the previously defined spatial variable **R** is replaced by **r** to signify that we are measuring the average dynamics in the plane of the cortex. Note that it is possible to extend the BOLD measurement to three dimensions to complement studies on laminar flows and cortical layer-specific activity [26]. However, this requires additional underlying theory to be developed and is beyond the scope of this study.

Transfer functions are vital in signal processing and control systems theory as they are powerful and convenient representations of the input–output response of a linear system [17]. Consider a transfer function *G*(*ω*) obtained through methods such as Fourier and Laplace transforms, which describes the response of a linear system to different frequencies *ω*. If *G*(*ω*) can be expressed as a ratio of two polynomials (i.e. the roots of *B*(*ω*) are poles of *G*(*ω*) and will give us an insight into the behaviour of the system for any input. Solving the linear dispersion relation corresponding to the poles will result in output solutions that represent the underlying modes of oscillation of the system. For example, if a pole of the form −*a* ± *bi* is found for *G*(*ω*), where *a* and *b* are positive constants, then we can immediately infer that one of the underlying modes of the system under study is a decaying sinusoid with decay rate *a* and oscillation frequency *b*. Other modes that can possibly arise from analysis of the poles are widely discussed in several control systems books such as [17]. Finally, the general response after an arbitrary excitation will be a weighted sum of all the modes.

In this section, we rewrite the linear BOLD transfer function in equation (2.10) as a ratio of two polynomials and calculate the poles. We decompose the resulting equation into components corresponding to each pole via the partial fraction method, thereby allowing us to derive the natural modes of the BOLD response.

Here we provide an analytic method to decompose the BOLD transfer function to pole components associated with its modes. We want to write *T _{Yζ}* as a linear sum of transfer functions

3.1

where *N* is the number of modes we can extract from *T _{Yζ}*. This enables us to characterize the separate contribution of the

3.2

such that the *j*th mode of *T _{Yζ}* follows the dispersion relation

3.3

and *a _{j}*(

We start by expressing equation (2.10) as a ratio of two polynomials *A*(*ω*) and *B*(*k*, *ω*) such that

3.4

where

3.5

and

3.6

We can write *A*(*ω*) in a compact form by defining real-valued constants

3.7

3.8

3.9

such that equation (3.5) becomes

3.10

On the other hand, *B*(*k*, *ω*) in equation (3.6) is a fifth-order polynomial, which implies that the BOLD transfer function has five frequency poles corresponding to five intrinsic modes (*N* = 5). To derive the exact form of the poles, we rewrite equation (3.6) as

3.11

where the dispersion relation in equation (3.3) is followed. Comparing equations (3.6) and (3.11), the complex frequency poles are

3.12

3.13

3.14

3.15

3.16

where *ω*_{3}, *ω*_{4} and *ω*_{5} are *k*-independent constants. We use equations (3.1), (3.2), (3.4), (3.10) and (3.11) to get the final expression for *a _{j}*(

3.17

where *δ _{jm}* is the Kronecker delta equal to 1 if

In §3.1, we analytically decomposed *T _{Yζ}* into a linear sum of component transfer functions

3.18

where

3.19

Many studies have shown that a linear section of the primary visual cortex will be stimulated if presented with a ring stimulus centred in the subject's visual field [13,27], owing to the natural retinotopic mapping between visual field and visual cortex [28]. Thus, a ring visual stimulus translates to a neural drive that only spatially varies perpendicular to a line (i.e. it is independent of the *y* direction) and results in a stronger BOLD response than its two-dimensional point-stimulus counterpart [13,15]. Hence, it is also worth calculating the BOLD response to a neural stimulus that is uniform in the *y*-dimension such that

3.20

where *k _{x}* and

3.21

We can then represent equation (3.21) as

3.22

where

3.23

is the *j*th component of the BOLD response in one dimension (i.e. to a line stimulus).

For a given neural drive *ζ*, the BOLD response *Y* depends on the properties of the transfer function *T _{Yζ}*(

The wave properties of the BOLD response *Y* depend on the intrinsic frequency characteristics of *T _{j}*(

Frequency characteristics of the component transfer functions for *v*_{β} = 1 mm s^{−1} and *Γ* = 1 s^{−1}. Each panel shows a contour plot in logarithmic scale of (*a*) |*T*_{1}(**k**, *ω*)|^{2}, (*b*) |*T*_{2}(**k**, *ω*)|^{2}, (*c*) |*T*_{3}(**k**, *ω* **...**

We observe the following relationships

4.1

and

4.2

signifying responses of the same magnitude that are mirror images of each other with respect to *ω* = 0. This remark can be further illustrated in coordinate space by taking the inverse Fourier transform where is the inverse Fourier transform operator, as plotted in figure 3. Figure 3 shows that where *F** denotes the complex conjugate of *F*, verifying that they represent travelling waves of the same speed but propagating in opposite directions. Also, we observe that *F*_{3}–*F*_{5} are local and non-propagating responses; with *F*_{3} and *F*_{4} having the same magnitude but opposite phases (because ), whereas *F*_{5} is real-valued having low amplitude.

Real and imaginary parts of *F*_{j}(**r**, *t*) for *v*_{β} = 1 mm s^{−1} and *Γ* = 1 s^{−1}. Each panel shows a contour plot normalized to the maximum of Re[*F*_{Yζ}(*r*, *t*)] as indicated by the colour bars. (*a*) Re[*F*_{1}(*r*, *t*)], (*b*) Re[*F*_{2}(*r*, *t* **...**

The broken lines in figure 2 show the maximums of |*T _{j}*(

4.3

where Re[·] denotes taking the real part and the explicit forms of *ω _{j}*(

Using parameters *v _{β}* = 1 mm s

In this section, we examine the spectral distribution of power in *T _{j}*(

4.4

and

4.5

In figure 4, we plot the temporal power as a function of frequency *f* = *ω*/2*π* for different pairs of *v _{β}* and

4.6

where Re[·] denotes taking the real part. For *j* = 3, 4, 5, where *ω _{j}* is independent of

4.7

which coincides with the results of our resonance calculation in equation (4.3). Upon substitution of nominal parameter values, we get
and which give the frequencies of the peaks of
and in figure 4. However, for *j* = 1, 2, where *ω _{j}* is

Temporal power on a logarithmic scale as a function of *f* = *ω*/2*π* for various combinations of *v*_{β} and *Γ*. Increasing colour shade from light to dark represents data for
and respectively. The parameters used in (*v*_{β} **...**

On the log–log scale used in the plot of spatial power versus spatial frequency *k* shown in figure 5, the general shape of for *Γ* = 1 s^{−1} is constant at low *k* and asymptotically approaches a power law of exponent –4 for higher *k* as expected from our theoretical formulae. The transition point between the flat and decreasing spectral regions moves to lower *k* as *v _{β}* increases, because the haemodynamic waves are able to propagate faster and are not limited to local dynamics. For

4.8

4.9

4.10

Upon substitution of nominal parameter values, equations (4.8)–(4.10) predict the frequencies of the peaks in figure 5. In §4.3, we study the mechanism behind the occurrence of maximum spatial power by relating to the resonances of *a _{j}*(

The form of *a _{j}*(

Functional forms of *ω*_{1}(*k*) and *ω*_{2}(*k*) versus *k* for various combinations of *v*_{β} and *Γ*. Each panel shows the real (solid lines) and imaginary (dashed lines) parts of *ω*_{1}(*k*) (black) and *ω*_{2}(*k*) (grey). The parameters **...**

We can easily calculate the value of *k*_{c} by setting *ω*_{1}(*k*_{c}) = *ω*_{2}(*k*_{c}). Substituting equations (3.12) and (3.13), the equality holds when

4.11

This result tells us that for *Γ* > *k _{z}v_{β}*,

Contour plot of *k*_{c} for various combinations of *v*_{β} and *Γ*. The solid line is a straight line with slope –*k*_{0} that divides the parameter space into real- and imaginary-valued *k*_{c}. Contours to the left side of the solid boundary line **...**

With the discovery of the boundary line in the *v _{β}* –

We now describe the behaviour of *a _{j}*(

Behaviour of *a*_{j}(*k*) as a function of *k* using *v*_{β} = 1 mm s^{−1} and different *Γ*. Decreasing colour shade from dark to light represents data for *a*_{1}, *a*_{2}, *a*_{3}, *a*_{4} and *a*_{5}, respectively. (*a*) |*a*_{j}| when *Γ* = 1.0 s^{−1}, (*b*) Arg[ **...**

On the other hand, for *v _{β}* = 1 mm s

In the previous sections, we analytically decomposed the total transfer function *T _{Yζ}*(

5.1

where *δ*(·) is the Dirac delta function. Once we have the Green function, we can immediately find the two-dimensional BOLD response *Y*(**r**, *t*) to an arbitrary stimulus *ζ*(**r**, *t*) by performing a convolution over space and time, with

5.2

where is the convolution operator and the integrals are over all space **r**′ and time *t*′. This is an alternative to the Fourier transform method of equation (2.11). The same principle is used to get the one-dimensional BOLD response *Y*(*x*, *t*) from the one-dimensional Green function *G*(*x*, *t*). In this section, we describe the integral forms of the Green function corresponding to each transfer function component. We also analytically derive closed forms of the Green function components for poles *ω*_{3}, *ω*_{4} and *ω*_{5}, which will significantly reduce the complexity of calculating the BOLD response in applications.

Here we describe how to calculate the Green function corresponding to *T _{j}*(

5.3

which shows that the Green function is just the inverse Fourier transform of the transfer function. Following the decomposition of *T _{Yζ}*(

5.4

where

5.5

Figure 9 shows the two-dimensional Green function components *G _{j}*(

Real and imaginary parts of the two-dimensional Green function components *G*_{j}**(****r**, *t*), where *r* = |**r**|, for *v*_{β} = 1 mm s^{−1} and *Γ* = 1 s^{−1}. Each panel shows a contour plot of the response normalized to the maximum of Re[*G*_{sum}(*r* **...**

Real and imaginary parts of the one-dimensional Green function components *G*_{j}(*x*, *t*) for *v*_{β} = 1 mm s^{−1} and *Γ* = 1 s^{−1}. Each panel shows a contour plot of the response normalized to the maximum of Re[*G*_{sum}(*x*, *t*)] as indicated **...**

Another interesting result is the non-zero response of *G*_{j}(**r**, *t*) at *t* = 0 and (figure 9*a–e*,*g*–*k*). These non-zero responses appear even though the neural stimulus used to get the Green functions is localized and non-zero only at *t* = 0 and *r* = 0 and seem to violate causality. However, we verify that they are not artefacts of the numerical implementation of equation (5.5). We then tried adding *G*_{1}(**r**, *t*)–*G*_{5}(**r**, *t*) and found that the responses at *t* = 0 and *r* ≠ 0 perfectly cancel each other to produce the expected zero response of *G*_{sum}(**r**, *t*) at *t* = 0 for any *r* ≠ 0 (figure 9*f*).

Finally, we emphasize that the behaviour of *G _{j}*(

Similarly, the Green function in one dimension *G*(*x*, *t*) is

5.6

where

5.7

Figure 10 shows the one-dimensional Green function components *G _{j}*(

The results we get from the two- and one-dimensional Green functions provide a comprehensive picture of the features we would expect from a BOLD response. We show both cases here to better appreciate that our theory can be compatible with future fMRI experiments, which may be measured in a preferred dimensionality (either two- or one dimension), as what we will illustrate in §6.

In §5.1, we derived the integral form of the Green function. Because we know the general form of the component transfer function *T _{j}*, we can reduce the complexity of equations (5.5) and (5.7) by analytically evaluating the

For the two-dimensional case, the transfer function *T _{j}*(

5.8

with *r* = |**r**| and *J*_{0}(*kr*) is the zeroth-order Bessel function of the first kind [29]. Upon substituting equation (3.2) into (5.8), we get

5.9

The *ω _{j}* poles in equations (3.12)–(3.16) lie in the lower half of the complex plane, so we evaluate the

5.10

where is the Heaviside function that expresses causality. Equation (5.10) is known as the Hankel transform and is less complex and faster to evaluate than the three-dimensional Fourier transform in equation (5.5). Many numerical implementations are available in the literature that can evaluate the Hankel transform efficiently, such as the method proposed in [30].

Path of integration to evaluate the complex integrals in: (*a*) equations (5.9) and (5.11), (*b*) equation (5.21) for the *x* > 0 case and (*c*) equation (5.21) for the *x* < 0 case. Each panel shows a semicircle contour in the complex plane to **...**

For the one-dimensional case, we substitute equation (3.2) into (5.7) and rearrange to get

5.11

As in the two-dimensional case, we evaluate the *ω* integral using the contour in figure 11*a* and simplify *G _{j}*(

5.12

Equation (5.12) is a single Fourier transform and is less complex and faster to evaluate than the two-dimensional Fourier transform in equation (5.7).

Equations (3.14)–(3.16) imply that for *j* = 3, 4, 5, the are independent of *k*, so that we can write This allows us to derive analytical closed forms for *G*_{3}–*G*_{5}. We begin by expressing equation (3.17) as a product of *k*-independent and *k*-dependent terms, as follows

5.13

where the term in the curly brackets is *k*-dependent. Substituting equations (3.12) and (3.13) into (5.13), we find

5.14

for *j* = 3, 4 and 5. If we define, the constant

5.15

we can rewrite equation (5.14) as

5.16

For the two-dimensional case, we substitute equation (5.16) into (5.10) and rearrange to obtain

5.17

Consider the instance when

5.18

where Re[·] denotes taking the real part. Equation (5.18) implies that An integral relation, only true for and is

5.19

where *K*_{0} is the zeroth-order modified Bessel function of the second kind [31]. Using equation (5.19), the closed-form solution of equation (5.17) becomes

5.20

For situations where condition in equation (5.18) is not satisfied, that is our derivation is still valid upon using the transformation

For the one-dimensional case, we substitute equation (5.16) into (5.12) and rearrange to obtain

5.21

For *x* > 0, we can evaluate equation (5.21) by using the Cauchy residue theorem and the contour shown in figure 11*b* to enclose the pole This results in

5.22

On the other hand, for *x* < 0, we use the contour in figure 11*c* to enclose the pole This results in

5.23

Combining the results of equations (5.22) and (5.23), the closed-form solution of equation (5.21) is

5.24

We emphasize that the analytical forms of equations (5.20) and (5.24) are only true for *j* = 3, 4 and 5, because they correspond to the *k*-independent poles, which is required for our derivation to work. Furthermore, we verified that for any *ν*_{β} and *Γ* value, equations (5.20) and (5.24) perfectly match the results when equations (5.5) and (5.7) are numerically integrated using fast Fourier transform, serving as a check of the accuracy of our results. We can use these analytical results to significantly reduce the complexity of determining the BOLD response to an arbitrary neural stimulus, because our remaining task is only to numerically solve for *G*_{1} and *G*_{2} instead of all five *G _{j}*.

We showed in §3.2 that decomposing the linear BOLD transfer function into components allowed us to recover the components *Y _{j}* of the linear BOLD response to a neural drive, each corresponding to a mode of frequency

In this section, we obtain the spatio-temporal haemodynamic response to a neural stimulus that is localized and low in amplitude. This kind of stimulus has been shown to produce a linear BOLD response on the human visual area V1. To test the model results experimentally, we compare our theory with fMRI data previously collected from a visual experiment [13]. In that study, the authors used a visual stimulus consisting of flickering concentric rings to measure the linear BOLD response on V1. They then used a wavefront phase estimation technique to quantitatively obtain the speed and damping rate of the haemodynamic waves observed for each subject. We use their estimations of *v _{β}* and

We model the stimulus mathematically as

6.1

where *σ _{x}* and

Experimental and simulated BOLD response to a spatio-temporal Gaussian stimulus. Each panel shows a contour plot in *x* (distance from centre of stimulus) and *t* (time from stimulus onset) space of the normalized response as indicated by the colour bars. **...**

We first note that the time to peak of the simulated BOLD response in figure 12*b* has a small offset relative to the experimental data in figure 12*a*. We attribute this temporal discrepancy to potential factors such as influences of astrocytic delays and other neurovascular coupling processes [32,33]. Incorporating these factors requires additional theoretical groundwork and is beyond the scope of this paper. However, we will implement model extensions related to astrocyte activity in a subsequent study to improve fits.

If we disregard the small temporal discrepancy, the general shape of the simulated response is similar to the result of the experiment especially in replicating the spatio-temporal width of the positive part of the signal. This demonstrates that the neural stimulus chosen is well grounded.

To study the components, we extended the spatial region of interest from [–5, 5] to [–10, 10] mm to be able to note obscured properties away from the origin. The responses show similar dynamics to the Green function components in figures 9 and and10,10, but with notable increases in spatial and temporal extent. In particular, *Y*_{1} and *Y*_{2} still correspond to travelling wave components, whereas *Y*_{3}–*Y*_{5} correspond to non-propagating components of the BOLD signal. Moreover, *Y*_{1}–*Y*_{4} have non-zero imaginary parts and follow the relations:
and where Re[·] and Im[·] denote taking the real and imaginary parts, respectively. These relations enable the imaginary parts of *Y*_{1}–*Y*_{4} to cancel each other producing real functions *Y*_{1} + *Y*_{2} and *Y*_{3} + *Y*_{4} resulting to a BOLD signal that is purely real. From this decomposition, we can clearly see that the local dynamics near *x* = 0 of the simulated BOLD response is dictated by the non-propagating components, whereas the spatial extent is highly influenced by the travelling wave component. This provides an exhaustive picture of the constituents of the BOLD response.

One of our goals is to demonstrate the potential to selectively enhance or suppress particular modes that underlie the BOLD response. We show in this section one way to enhance the relative amplitude of a BOLD response component by using a moving spatio-temporal sinusoid as neural stimulus. We mathematically model the neural stimulus as

6.2

where *k*_{s} and *ω*_{s} are the real-valued spatial and temporal frequencies of the stimulus, respectively. In frequency space, the Fourier transform of equation (6.2) is a combination of Dirac delta functions

6.3

In order to get the *j*th component of the BOLD response, we either convolve equation (6.2) with *G _{j}*(

6.4

We further simplify equation (6.4) by substituting equation (3.2) and exploiting the fact that *a _{j}*(

6.5

where *a _{j}*(

6.6

and we can separate into its real and imaginary parts, with

6.7

We can see from equation (6.6) that there are two ways to enhance the response of *Y _{j}*(

We can approximate the degree of enhancement of the second method by taking the ratio of the absolute value of equation (6.6) when to its value when the stimulus is almost stationary in time; *ω*_{s} is very small or An approximate value of this ratio at positions and times where and is

6.8

which always greater than 1 for any *j*. Note that *ratio* is a constant for *j* = 3, 4 and 5, because *ω*_{3}, *ω*_{4} and *ω*_{5} are constants. For *j* = 1 and 2, for very weak or very strong damping of *ω*_{1}(*k*_{s}) and *ω*_{2}(*k*_{s}). This proves that a sinusoidal stimulus is able to enhance the response of the *j*th component of the BOLD response, as long as its frequency is tuned appropriately. This illustrative application will be implemented in future works to test its accuracy.

We have analysed the transfer function from a physiologically based poroelastic haemodynamic model [13–15], to provide a deeper understanding of the linear BOLD response to neural activity. In particular, we have developed quantitative tools to uncover the fundamental components that contribute to the overall structure of the haemodynamic response to a localized impulse input. The main results of this study are as follows:

- (i)We analytically decomposed the linear BOLD transfer function into components corresponding to frequency poles via partial fractions. The poles are responsible for resonances of the transfer function and are associated with modes of the BOLD response. Each mode has a distinct dispersion relation, as shown in equation (3.3).
- (ii)The derived component transfer functions were characterized in terms of frequency responses, power spectra and calculation of resonances. We found from the frequency responses that the BOLD transfer function consists of three components that are non-propagating and two components that represent waves travelling in opposite directions. Each component has resonance frequencies that we calculated analytically in equation (4.3) and illustrated in figure 2. Analysis of spatial and temporal power spectra show that the damping rate
*Γ*restricts significant responses in*r*–*t*space to a localized region because it increases flow losses owing to blood drainage. We also found that the components have peak values tuned to selected frequencies and ; exact values of which are shown in equations (4.6)–(4.10). The origin of the peak values is related to the properties of*a*(_{j}*k*). Because*a*(_{j}*k*) is proportional to the inverse of the differences of*ω*(_{j}*k*) pairs, critical spatial frequencies*k*_{c}where the equality*ω*_{1}(*k*) =*ω*_{2}(*k*) is achieved can lead to resonances that enhance the responses. We further showed how*k*_{c}changes, from real to imaginary values, as a function of propagation speed*v*and damping rate_{β}*Γ*. - (iii)We used the transfer function components to derive and decompose the BOLD response to a localized neural impulse—i.e. the Green function or stHRF. Our technique analytically decomposes the Green function into component responses that can be related to different physiological phenomena. This provides an efficient method to calculate the linear BOLD response to any kind of neural activity that is less complex and faster to implement than direct Fourier transform methods.
- (iv)The analytical results were used to demonstrate how the BOLD response to a spatio-temporal Gaussian neural stimulus can be decomposed into response components. For spatial and temporal widths of the stimulus matching the analysis of an experiment by [13], we decomposed the linear BOLD response into travelling waves and non-propagating components, which enabled us to clearly identify the origin of the shape and characteristics of the BOLD response in terms of underlying physiology.
- (v)An illustrative application was provided to show how a moving sinusoidal neural stimulus can be used to alter the effect of different components of the BOLD response. This was done by tuning the spatial and temporal frequencies of the stimulus to the preferred resonance frequencies of particular components. We estimated the potential increase in the response and quantitatively showed that our method is capable of enhancing the component corresponding to the chosen resonance frequency.

Overall, our present work establishes new methods to quantitatively understand the spatio-temporal properties of haemodynamic responses to various neural stimuli. It provides basis for further analyses and predictions for comparison with experiment. For example, it has the potential to aid in designing experiment protocols that can selectively enhance or suppress spatio-temporal resonances to increase or decrease signal amplitude of some responses. One could potentially design a visual stimulus, such as a contrast-reversing grating with sinusoidal contrast modulation [34], to elicit a neural response similar to the example in §6.2 and highlight the corresponding BOLD response component. In addition, because we found that the BOLD response consists of wave-like and non-propagating components, our results can also be used to design experiments to improve measurement of vascular properties such as *v _{β}* and

Another promising application is to derive the modes underlying the BOLD response to different types of stimulus. A recent study by [35] developed a method, based on the model we used, to deconvolve fMRI data and extract spatio-temporal patterns of neural activity. The formulated deconvolution technique would be able to separate neural and haemodynamic contributions on fMRI maps to gain insights of various physiological phenomena such as detection of stimulus orientation in human primary visual cortex [36]. This is crucial to obtain accurate neural activity maps that can disambiguate real neural dynamics from spatially distributed connectivity delays [37]. As advances in fMRI scanner hardware continue to improve, we can use the deconvolution technique to obtain an estimate of neural activity *ζ*. The estimated *ζ* can then be convolved with the Green function components *G _{j}* to calculate the response modes

The authors thank Thomas Lacy for useful discussions.

P.A.R. proposed the decomposition method and showed substantial feasibility. J.C.P. carried out most subsequent theoretical work and all numerical work. All authors performed analysis and interpretation of results and drafted the manuscript. All authors gave final approval for publication.

The authors declare no competing interests.

This work was supported by the Australian Research Council Center of Excellence for Integrative Brain Function (ARC Center of Excellence Grant CE140100007), the Australian Research Council Laureate Fellowship Grant (FL140100025), and the Australian Research Council Discovery Project (DP130100437).

1. Kwong KK, et al.
1992.
Dynamic magnetic resonance imaging of human brain activity during primary sensory stimulation. Proc. Natl Acad. Sci. USA. 89, 5675–5679. (doi:10.1073/pnas.89.12.5675) [PubMed]

2. Ogawa S, Tank DW, Menon R, Ellermann JM, Kim SG, Merkle H, Ugurbil K
1992.
Intrinsic signal changes accompanying sensory stimulation: functional brain mapping with magnetic resonance imaging. Proc. Natl Acad. Sci. USA
89, 5951–5955. (doi:10.1073/pnas.89.13.5951) [PubMed]

3. Hennig J, Speck O, Koch MA, Weiller C
2003.
Functional magnetic resonance imaging: a review of methodological aspects and clinical applications. J. Magn. Reson. Imaging
18, 1–15. (doi:10.1002/jmri.10330) [PubMed]

4. Buxton RB.
2009.
Introduction to functional magnetic resonance imaging: principles and techniques. Cambridge, UK: Cambridge University Press.

5. Friston KJ.
2009.
Modalities, modes, and models in functional neuroimaging. Science
326, 399–403. (doi:10.1126/science.1174521) [PubMed]

6. Buxton RB, Wong EC, Frank LR
1998.
Dynamics of blood flow and oxygenation changes during brain activation: the balloon model. Magn. Reson. Med. 39, 855–864. (doi:10.1002/mrm.1910390602) [PubMed]

7. Friston KJ, Mechelli A, Turner R, Price CJ
2000.
Nonlinear responses in fMRI: the balloon model, Volterra kernels, and other hemodynamics. Neuroimage
12, 466–477. (doi:10.1006/nimg.2000.0630) [PubMed]

8. Buxton RB, Uludauğ K, Dubowitz DJ, Liu TT
2004.
Modeling the hemodynamic response to brain activation. Neuroimage
23, S220–S233. (doi:10.1016/j.neuroimage.2004.07.013) [PubMed]

9. 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. (doi:10.1016/j.neuroimage.2009.09.059) [PMC free article] [PubMed]

10. Engel SA, Glover GH, Wandell BA
1997.
Retinotopic organization in human visual cortex and the spatial precision of functional MRI. Cereb. Cortex
7, 181–192. (doi:10.1093/cercor/7.2.181) [PubMed]

11. Schmuel A, Yacoub E, Chaimow D, Logothetis NK, Ugurbil K
2007.
Spatio-temporal point-spread function of functional MRI signal in human gray matter. Neuroimage
35, 539–552. (doi:10.1016/j.neuroimage.2006.12.030) [PMC free article] [PubMed]

12. Friston KJ.
1995.
Regulation of rCBF by diffusible signals: an analysis of constraints on diffusion and elimination. Hum. Brain Mapp. 3, 56–65. (doi:10.1002/hbm.460030106)

13. Aquino KM, Schira MM, Robinson PA, Drysdale PM, Breakspear MJ
2012.
Hemodynamic traveling waves in human visual cortex. PLoS Comput. Biol. 8, e1002435 (doi:10.1371/journal.pcbi.1002435) [PMC free article] [PubMed]

14. Drysdale PM, Huber JP, Robinson PA, Aquino KM
2010.
Spatiotemporal BOLD dynamics from a porous elastic hemodynamic model. J. Theor. Biol. 265, 524–534. (doi:10.1016/j.jtbi.2010.05.026) [PubMed]

15. Aquino KM, Robinson PA, Drysdale PM
2014.
Spatiotemporal hemodynamic response functions derived from physiology. J. Theor. Biol. 347, 118–136. (doi:10.1016/j.jtbi.2013.12.027) [PubMed]

16. Wang HF.
2000.
Theory of linear poroelasticity with applications to geomechanics and hydrogeology, 1st edn
Princeton, NJ: Princeton University Press.

17. Shinners SM.
1998.
Modern control system theory and design, 2nd edn
New York, NY: John Wiley & Sons.

18. Mansouri R, Bettayeb M, Djennoune S
2010.
Approximation of high order integer systems by fractional order reduced-parameters models. Math. Comp. Model. 51, 53–62. (doi:10.1016/j.mcm.2009.07.018)

19. Grubb RL, Raichle ME, Eichling JO, Ter-Pogossian MM
1974.
The effects of changes in PaCO_{2} on cerebral blood volume, blood flow, and vascular mean transit time. Stroke
5, 630–639. (doi:10.1161/01.STR.5.5.630) [PubMed]

20. Stephan KE, Weiskopf N, Drysdale PM, Robinson PA, Friston KJ
2007.
Comparing hemodynamic models with DCM. Neuroimage
38, 387–401. (doi:10.1016/j.neuroimage.2007.07.040) [PMC free article] [PubMed]

21. Robinson PA, Drysdale PM, der Merwe HV, Kyriakou E, Rigozzi MK, Germanoska B, Rennie CJ
2006.
BOLD responses to stimuli: dependence on frequency, stimulus form, amplitude, and repetition rate. Neuroimage
31, 585–599. (doi:10.1016/j.neuroimage.2005.12.026) [PubMed]

22. Friston KJ, Harrison L, Penny W
2003.
Dynamic causal modelling. Neuroimage
19, 1273–1302. (doi:10.1016/S1053-8119(03)00202-7) [PubMed]

23. Trudnowski RJ, Rico RC
1974.
Specific gravity of blood and plasma at 4 and 37°C. Clin. Chem. 20, 615–616. [PubMed]

24. Huppert TJ, Hoge RD, Diamond SG, Franceschini MA, Boas DA
2006.
A temporal comparison of BOLD, ASL, and NIRS hemodynamic responses to motor stimuli in adult humans. Neuroimage
29, 368–382. (doi:10.1016/j.neuroimage.2005.08.065) [PMC free article] [PubMed]

25. Duvernoy HM, Delon S, Vannson JL
1981.
Cortical blood vessels of the human brain. Brain Res Bull. 7, 519–579. (doi:10.1016/0361-9230(81)90007-1) [PubMed]

26. Olman CA, Harel N, Feinberg DA, He S, Zhang P, Ugurbil K, Yacoub E et al.
2012.
Layer-specific fMRI reflects different neuronal computations at different depths in human V1. PLoS ONE
7, e32536 (doi:10.1371/journal.pone.0032536) [PMC free article] [PubMed]

27. Tootell RB, Silverman MS, Switkes E, Valois RLD
1982.
Deoxyglucose analysis of retinotopic organization in primary striate cortex. Science
218, 902–904. (doi:10.1126/science.7134981) [PubMed]

28. Schira MM, Tyler CW, Breakspear MJ, Spehar B
2009.
The foveal confluence in human visual cortex. J. Neurosci. 29, 9050–9058. (doi:10.1523/JNEUROSCI.1760-09.2009) [PubMed]

29. Gradshteyn IS, Ryzhik IM, Jeffrey A, Zwillinger D
2007.
Table of integrals, series, and products, 7th edn
New York, NY: Academic Press.

30. Guizar-Sicairos M, Gutiérrez-Vega JC
2004.
Computation of quasi-discrete Hankel transforms of integer order for propagating optical wave fields. J. Opt. Soc. Am. A
21, 53–58. (doi:10.1364/JOSAA.21.000053) [PubMed]

31. Prudnikov AP, Brychkov YA, Marichev OI
1986.
Integrals and series volume 2: special functions. Amsterdam, The Netherlands: Gordon and Breach.

32. DeYoe EA, Bandettini P, Neitz J, Miller D, Winans P
1994.
Functional magnetic resonance imaging (FMRI) of the human brain. J. Neurosci. Methods
54, 171–187. (doi:10.1016/0165-0270(94)90191-0) [PubMed]

33. Masamoto K, Unekawa M, Watanabe R, Toriumi H, Takuwa H, Kawaguchi H et al.
2015.
Unveiling astrocytic control of cerebral blood flow with optogenetics. Sci. Rep. 5, 11455 (doi:10.1038/srep11455) [PMC free article] [PubMed]

34. Ferree TC, Wade AR, Calvisi ML, Szeri AJ, Liley DTJ, Inglis BA et al.
2007.
Experimental study of EEG and BOLD responses to sinusoidal contrast modulation. In *Engineering in Medicine and Biology Workshop*, *2007 IEEE Dallas*, *TX, 11–12 November*, pp. 12–15. Piscataway, NJ: IEEE.

35. Aquino KM, Robinson PA, Schira MM, Breakspear MJ
2014.
Deconvolution of neural dynamics from fMRI data using a spatiotemporal hemodynamic response function. Neuroimage
94, 203–215. (doi:10.1016/j.neuroimage.2014.03.001) [PubMed]

36. Haynes JD, Rees G
2005.
Predicting the orientation of invisible stimuli from activity in human primary visual cortex. Nat. Neurosci. 8, 686–691. (doi:10.1038/nn1445) [PubMed]

37. Chang C, Thomason ME, Glover GH
2008.
Mapping and correction of vascular hemodynamic latency in the BOLD signal. Neuroimage
43, 90–102. (doi:10.1016/j.neuroimage.2008.06.030) [PMC free article] [PubMed]

Articles from Journal of the Royal Society Interface are provided here courtesy of **The Royal Society**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |