## Abstract

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.

## 1. Introduction

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.

## 2. Spatio-temporal haemodynamic model

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].

### 2.1. Dynamical equations

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.1where *κ*, *γ* and *F*_{0} are the blood flow signal decay rate, the flow-dependent elimination constant and the resting blood inflow, respectively.

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.2where *D* quantifies damping owing to effective blood viscosity, *ρ _{f}* is the blood density,

*c*

_{1}is the pressure coupling constant,

*δ*(

*z*) is the Dirac delta function expressing the inflows or outflows at the cortical surface

*z*= 0,

*c*is the blood outflow constant, and

_{P}*P*(

**R**,

*t*) is the average pressure. Equation (2.2) predicts propagation of Ξ waves with serving as the spatial coupling between adjacent sites in the tissue owing to pore pressure gradients. In addition,

*P*(

**R**,

*t*) is related to Ξ via the constitutive equation 2.3where

*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.4where **v**(**R**, *t*) is the average blood velocity, **v*** _{F}*(

**R**,

*t*) is the inflow blood velocity,

**v**

*(*

_{P}**R**,

*t*) is the outflow blood velocity,

*ψ*is the ratio of haemoglobin concentration to blood density and

*η*is the fractional rate of oxygen consumption. The first term on the right-hand side relates the rate of change of dHb to the amount of flux going to adjacent sites, the second term is the conversion of dHb to oxygenated haemoglobin (

*ψ*Ξ −

*Q*) at a rate

*η*, and the last term represents the rate of reduction of dHb owing to blood outflow at the boundary. On the other hand,

**v**

*(*

_{F}**R**,

*t*) and

**v**

*(*

_{P}**R**,

*t*) are derived from the boundary conditions imposed on

*Q*(

**R**,

*t*) and follow the relations 2.5and 2.6whereas 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.8where *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*).

### 2.2. Linear blood oxygen-level dependent transfer function

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

*θ*(i.e.

*F*, Ξ,

*Q*,

*Y*) per unit change in neural drive

*ζ*at the same

**k**and

*ω*. As an example, for

*θ*=

*Y*, the transfer function

*T*that relates the BOLD signal to the neural drive is 2.9where

_{Yζ}*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.10where is the spatial wavenumber in the plane of the cortex, *C _{z}* is the outflow normalization constant,

*k*is the effective spatial frequency,

_{z}*ω*

_{f}is the natural frequency of flow response,

*τ*is the haemodynamic transit time and

*v*and

_{β}*Γ*are the free parameters corresponding to the wave propagation speed and damping rate, respectively [13,15]. All constants are either obtained from known physiological estimates or derived from steady-state behaviour [15]; their nominal values from physiology are summarized in table 1.

### 2.3. Calculation of the blood oxygen-level dependent response

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.11where 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.

## 3. Natural modes of the blood oxygen-level dependent response

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.

### 3.1. Pole decomposition of the blood oxygen-level dependent transfer function

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

*T*(

_{j}**k**,

*ω*), i.e. 3.1where

*N*is the number of modes we can extract from

*T*. This enables us to characterize the separate contribution of the

_{Yζ}*j*th mode to the overall BOLD response. In addition, we want

*T*to be related to a temporal frequency pole

_{j}*ω*(

_{j}*k*) with general form 3.2such that the

*j*th mode of

*T*follows the dispersion relation 3.3and

_{Yζ}*a*(

_{j}*k*) is a

*k*-dependent quantity that needs to be determined. Equations (3.1) and (3.2), when combined, decompose the BOLD transfer function via the partial fraction method. The rest of the section is devoted to the derivation of the number of modes

*N*, the poles

*ω*(

_{j}*k*), and the quantities

*a*(

_{j}*k*).

We start by expressing equation (2.10) as a ratio of two polynomials *A*(*ω*) and *B*(*k*, *ω*) such that
3.4where
3.5and
3.6We can write *A*(*ω*) in a compact form by defining real-valued constants
3.7
3.8
3.9such that equation (3.5) becomes
3.10On 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.11where 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.16where *ω*_{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}*(

*k*) as 3.17where

*δ*is the Kronecker delta equal to 1 if

_{jm}*j*=

*m*and 0 otherwise. We see from the form of

*a*(

_{j}*k*) that each

*T*depends on all

_{j}*m*≠

*j*poles, implying that the modes of the BOLD response will have inter-related properties.

### 3.2. Decomposition of the blood oxygen-level dependent response

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

*T*. This result when substituted in equation (2.11) leads to the decomposition of two-dimensional BOLD response

_{j}*Y*(

**r**,

*t*) as a linear sum of five response components

*Y*(

_{j}**r**,

*t*). Mathematically, this is demonstrated as 3.18where 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.20where *k _{x}* and

*k*are the spatial frequencies in the

_{y}*x-*and

*y*-directions, respectively. Substituting equation (3.20) into (2.11) will allow us to drop the

*y*argument and

*Y*will only depend on the one-dimensional distance

*x*perpendicular to the centre of stimulus given by 3.21

We can then represent equation (3.21) as
3.22where
3.23is the *j*th component of the BOLD response in one dimension (i.e. to a line stimulus).

## 4. Features of components *T*_{j}(**k**,*ω*)

_{j}

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

**k**,

*ω*). The decomposition derived in §3 allows us to isolate the contribution of each component

*T*(

_{j}**k**,

*ω*) to the overall characteristics of

*T*(

_{Yζ}**k**,

*ω*), providing a detailed understanding of the BOLD response. It also aids in formulating possible methods to selectively enhance or suppress specific modes of the BOLD response depending on which pole

*ω*(

_{j}*k*) we want to focus on. In this section, we describe the frequency characteristics of

*T*(

_{j}**k**,

*ω*) and calculate its corresponding temporal and spatial power spectra. We analytically derive the frequencies of resonances in the power spectra and associate them with the singularities of

*a*(

_{j}*k*).

### 4.1. Spatio-temporal frequency properties

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

**k**,

*ω*). We show in figure 2 the magnitude of

*T*using

_{j}*v*= 1 mm s

_{β}^{−1}and

*Γ*= 1 s

^{−1}for different spatio-temporal frequencies

*k*and

*f*=

*ω*/2

*π*.

We observe the following relationships
4.1and
4.2signifying 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.

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

**k**,

*ω*)|

^{2}representing resonances that indicate where the BOLD response

*Y*will preferentially oscillate. For each

*k*value, the broken lines predict a corresponding temporal resonance frequency

*f*

_{j}_{,res}=

*ω*

_{j}_{,res}/2

*π*. The general form or value of

*ω*

_{j}_{,res}is exactly solved by taking the partial derivative of |

*T*(

_{j}**k**,

*ω*)|

^{2}with respect to

*ω*and equating it to zero. This results in 4.3where Re[·] denotes taking the real part and the explicit forms of

*ω*(

_{j}*k*) are given in equations (3.12)–(3.16).

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

^{−1}and

*Γ*= 1 s

^{−1}, and nominal values for the constants, the relation is valid, such that and For

*k*-values where

*ω*

_{1,res}and

*ω*

_{2,res}vary almost linearly with

*k*as shown by the broken lines in figure 2

*a*,

*b*. In addition, using

*ω*= 0.56 s

_{f}^{−1}(table 1), the resonances

*ω*

_{3,res}/2

*π*= −0.089,

*ω*

_{4,res}/2

*π*= +0.089 and

*ω*

_{5,res}= 0 Hz are the specific frequencies of the broken lines in figure 2

*c*,

*d*,

*e*, respectively. Moreover, the resonances for the total transfer function in figure 2

*f*occur at all resonance frequencies observed for its components.

### 4.2. Power spectra

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

**k**,

*ω*). We do this by calculating the temporal power spectrum and spatial power spectrum as follows 4.4and 4.5In figure 4, we plot the temporal power as a function of frequency

*f*=

*ω*/2

*π*for different pairs of

*v*and

_{β}*Γ*. We observe that each

*T*(

_{j}**k**,

*ω*) has high temporal power at low frequencies, characteristic of a low-pass filter and consistent with the findings of experiments [14] and the blood volume fluctuation transfer function predicted by the balloon model [21]. We also see that by independently increasing

*v*and

_{β}*Γ*, the amplitudes of all spectra decreases. However, higher

*Γ*broadens and but does not change the shape of and Moreover, it increases the low-frequency (near

*f*= 0 Hz) value of because it causes the response to be more localized, thus it has a higher amplitude at its centre. Furthermore, each reaches peak values at certain resonance frequencies The value of can generally be calculated by taking the derivative of with respect to

*ω*and equating it to zero, which results in 4.6where Re[·] denotes taking the real part. For

*j*= 3, 4, 5, where

*ω*is independent of

_{j}*k*, the solution for equation (4.6) can be obtained analytically such that 4.7which 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

*ω*is

_{j}*k*-dependent, equation (4.6) is not analytically tractable and we need to use numerical methods to get the exact values of and However, for cases when such as for

*v*–

_{β}*Γ*pairs (1.0 mm s

^{–1}, 0.6 s

^{–1}) and (2.5 mm s

^{–1}, 0.6 s

^{–1}), which is evident in figure 4

*d*,

*e*. We further note that and move to higher values as

*v*and

_{β}*Γ*are increased.

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

*Γ*= 0.6 s

^{−1}combined with

*v*= 1 mm s

_{β}^{−1}or 2.5 mm s

^{–1}, the constant-to-power law transition still occurs for and However, for and we see a constant spectrum for lower

*k*followed by a maximum at an intermediate

*k*before it decreases with an asymptotic power law of exponent –4 for higher

*k*. The frequency of the maximum, if it exists, can generally be calculated by taking the derivative of with respect to

*k*and equating it to zero, which yields 4.8 4.9 4.10Upon 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}*k*).

### 4.3. Singularity of *a*_{j}(*k*)

_{j}

The form of *a _{j}*(

*k*) in equation (3.17) is proportional to [

*ω*(

_{j}*k*) −

*ω*(

_{m}*k*)]

^{−1}which will produce singularities if

*ω*(

_{j}*k*) =

*ω*(

_{m}*k*) for

*j*≠

*m*. This is not applicable for

*j*,

*m*= 3, 4, 5, because

*ω*

_{3},

*ω*

_{4}and

*ω*

_{5}are

*k*-independent constants having different values. However, for the case of

*j*,

*m*= 1, 2,

*ω*

_{1}(

*k*) and

*ω*

_{2}(

*k*) are

*k*-dependent quantities that have the same value at a critical spatial frequency

*k*

_{c}. In figure 6

*d,e*, we simultaneously plot the real and imaginary parts of

*ω*

_{1}(

*k*) and

*ω*

_{2}(

*k*) and show that a solution for

*k*

_{c}exists when

*Γ*= 0.6 s

^{−1}and

*v*≤ 2.5 mm s

_{β}^{−1}. We also observe that the real parts of

*ω*

_{1}(

*k*) and

*ω*

_{2}(

*k*) remain zero for

*k*≤

*k*

_{c}, showing that the poles are purely imaginary and correspond to damped responses. However, in figures 6

*a*,

*b,c,f*, where

*v*and

_{β}*Γ*are varied, for all

*k*. These results imply that the existence of

*k*

_{c}depends strongly on the parameters

*v*and

_{β}*Γ*.

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.11This result tells us that for *Γ* > *k _{z}v_{β}*,

*k*

_{c}is real such that

*ω*

_{1}(

*k*) can have the same value with

*ω*

_{2}(

*k*) in the real

*k*-space; thus, singularity of

*a*(

_{j}*k*) can be observed. To obtain an idea of the range of values

*k*

_{c}may have, we show in figure 7 the values of

*k*

_{c}for various combinations of

*v*and

_{β}*Γ*. Figure 7 shows a solid boundary line that separates the

*v*–

_{β}*Γ*parameter space into two regions where

*k*

_{c}is either real or imaginary. The equation of the line is which can be obtained by imposing

*k*

_{c}= 0 and substituting the nominal values from table 1. The importance of this line is that we can immediately discern the nature of

*k*

_{c}for a given pair of

*v*and

_{β}*Γ*even without calculating its explicit value. If the pair falls to the left side of the boundary line,

*k*

_{c}is immediately real, and we would expect singularities in

*a*(

_{j}*k*).

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

*Γ*parameter space, it is easy to see that

*Γ*= 1.0 s

^{−1}combined with any

*v*falls to the right side of the boundary line, thus leading to imaginary-valued

_{β}*k*

_{c}. This will result in non-intersecting

*ω*

_{1}(

*k*) and

*ω*

_{2}(

*k*), as shown in figures 6

*b*,

*c*. The same analysis is true when

*Γ*= 0.6 s

^{−1}is paired with

*v*= 5 mm s

_{β}^{−1}, as verified in figure 6

*f*. However, when

*Γ*= 0.6 s

^{−1}is combined with any

*v*< 2.58 mm s

_{β}^{−1}, they fall to the left side of the boundary line, thus

*k*

_{c}is real and

*ω*

_{1}(

*k*

_{c}) =

*ω*

_{2}(

*k*

_{c}), as shown in figure 6

*d*,

*e*.

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

*k*) in figure 8 by calculating its modulus and argument, using

*v*= 1 mm s

_{β}^{−1}and different damping rates. The damping rates chosen are

*Γ*= 0.6 and 1.0 s

^{–1}to investigate the properties of

*a*(

_{j}*k*) when

*k*

_{c}is either real or imaginary. In figure 8

*a*,

*b*, where

*v*= 1 mm s

_{β}^{−1}and

*Γ*= 1.0 s

^{−1},

*k*

_{c}is imaginary, thereby leading to monotonic changes in |

*a*(

_{j}*k*)|. Also, we see that |

*a*

_{1}(

*k*)| overlaps with |

*a*

_{2}(

*k*)|, as well as |

*a*

_{3}(

*k*)| with |

*a*

_{4}(

*k*)|, which can be attributed to the similarity in the forms of their corresponding poles. The discontinuous jumps in the plot of Arg[

*a*] are artefacts of the restriction of the range of Arg to [−

_{j}*π*,

*π*), and can easily be removed by translating the range to [0, 2

*π*).

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

^{−1}combined with

*Γ*= 0.6 s

^{−1}(figure 8

*c*,

*d*),

*k*

_{c}is real. This is why we see a sharp resonance for |

*a*

_{1}(

*k*)| and |

*a*

_{2}(

*k*)| at where is defined in equation (4.8). We also observe a resonance for |

*a*

_{3}(

*k*)| and |

*a*

_{4}(

*k*)| at with from equation (4.9). These resonances are responsible for the occurrence of the peaks in the spatial power spectra in figure 5, because is proportional to

*a*(

_{j}*k*).

## 5. Blood oxygen-level dependent response to a localized neural activity: Green function

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

**k**,

*ω*) into components and analysed their characteristics, in order to retrieve and describe the underlying modes of oscillation of the BOLD response

*Y*. We now want to explore the effects of the nature of these modes on the haemodynamic response. However, instead of calculating the effect of

*T*(

_{j}**k**,

*ω*) on the BOLD response using the Fourier transform method of §2.3, we relate it to the Green function, also called the stHRF. In general, the Green function characterizes the response of a system to the presence of an impulse source. Because any distribution of source can be written as a combination of impulse sources, knowing how the system reacts to an impulse source provides significant information as to how the system will react to a distribution of source. In this paper's context, the two-dimensional Green function

*G*(

**r**,

*t*) is the BOLD response to neural activity that is localized in space and time, i.e. 5.1where

*δ*(·) 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.2where ⊗ 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.

### 5.1. General form of the Green function

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

**k**,

*ω*). We start by considering that the Fourier transform of equation (5.1) is

*ζ*(

**k**,

*ω*) = 1. When substituted in equation (2.11), this gives the general integral form of the two-dimensional Green function

*G*(

**r**,

*t*) as 5.3which shows that the Green function is just the inverse Fourier transform of the transfer function. Following the decomposition of

*T*(

_{Yζ}**k**,

*ω*) and

*Y*(

**r**,

*t*) in equations (3.1) and (3.18), respectively, we can express the Green function as the linear sum 5.4where 5.5

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

**r**,

*t*) in coordinate space using and We first note that

*G*

_{1}(

**r**,

*t*) and

*G*

_{2}(

**r**,

*t*) have wave-like responses propagating in a small spatial range with speeds comparable to

*v*. The structure of the wavefronts will be seen more clearly in one dimension (figure 10

_{β}*a,b*) as the waves are able to propagate further in space because they only spread in one direction, as opposed to two for two dimensions. Thus, the effective amplitude of a one-dimensional wave will fall off proportional to instead of for two dimensions. On the other hand,

*G*

_{3}(

**r**,

*t*),

*G*

_{4}(

**r**,

*t*) and

*G*

_{5}(

**r**,

*t*) remain as non-propagating responses, agreeing with the results of §4.1. Surprisingly,

*G*

_{1}(

**r**,

*t*)–

*G*

_{4}(

**r**,

*t*) have non-zero imaginary parts, whereas

*G*

_{5}(

**r**,

*t*) is negative and purely real. This does not follow a simple expectation that the components must be real-valued functions, because the total Green function

*G*

_{sum}(

**r**,

*t*) is real. However, upon closer inspection, we find that the components follow the relations: and where Re[·] and Im[·] denote taking the real and imaginary parts, respectively. Thus, the combined responses and produce real-valued functions and can potentially be related to observable haemodynamic quantities. These results account for why the total BOLD response

*G*

_{sum}(

**r**,

*t*) remains real (figures 9

*f,l*).

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}*(

**r**,

*t*) we see in figure 9 is general but the exact spatio-temporal variation depends strongly on the chosen properties of the vasculature, particularly the wave propagation speed

*v*and damping rate

_{β}*Γ*; the spatial range of the responses increases for higher

*v*and decreases for higher

_{β}*Γ*.

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

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

*x*,

*t*) in coordinate space using

*v*= 1 mm s

_{β}^{−1}and

*Γ*= 1 s

^{−1}. We see that

*G*(

_{j}*x*,

*t*) has a symmetric activity on both sides of

*x*= 0 but the shape is generally similar to its two-dimensional counterpart

*G*(

_{j}**r**,

*t*). However, the responses are stronger and have larger spatio-temporal extent, in agreement with the findings of [13,15]. We also note that

*G*

_{1}(

*x*,

*t*) and

*G*

_{2}(

*x*,

*t*) still correspond to travelling waves with propagation speeds comparable to

*v*, because the wavefronts move almost parallel to the guide line shown in figure 10

_{β}*a*,

*b*, whereas

*G*

_{3}(

*x*,

*t*),

*G*

_{4}(

*x*,

*t*) and

*G*

_{5}(

*x*,

*t*) are non-propagating responses, similar to the findings of figure 9. More importantly, we still see the unique characteristics of the components such as (i) non-zero imaginary parts for

*G*

_{1}(

*x*,

*t*)–

*G*

_{4}(

*x*,

*t*) and (ii) the paradoxical non-zero response of

*G*(

_{j}*x*,

*t*) at

*t*= 0 and

*x*≠ 0. The imaginary parts then cancel each other such that

*G*

_{1}(

*x*,

*t*) +

*G*

_{2}(

*x*,

*t*) and

*G*

_{3}(

*x*,

*t*) +

*G*

_{4}(

*x*,

*t*) are real. In addition, combining all

*G*(

_{j}*x*,

*t*), the response at

*t*= 0 and

*x*≠ 0 vanishes as expected for the BOLD response (figure 10

*f*,

*l*).

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.

### 5.2. Simplifying the integral form of the Green functions

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

*ω*integral. In this section, we show how this can be done using contour integration.

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

**k**,

*ω*) is radially symmetric in

**k**, thus

*G*(

_{j}**r**,

*t*) will also be radially symmetric such that equation (5.5) becomes 5.8with

*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

*ω*integral using Cauchy residue theorem and the contour shown in figure 11

*a*to enclose

*ω*(

_{j}*k*). This simplifies

*G*(

_{j}*r*,

*t*) to 5.10where 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].

For the one-dimensional case, we substitute equation (3.2) into (5.7) and rearrange to get
5.11As in the two-dimensional case, we evaluate the *ω* integral using the contour in figure 11*a* and simplify *G _{j}*(

*x*,

*t*) to 5.12Equation (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).

### 5.3. Analytic derivation of the Green function for *j* = 3, 4 and 5

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.13where the term in the curly brackets is *k*-dependent. Substituting equations (3.12) and (3.13) into (5.13), we find
5.14for *j* = 3, 4 and 5. If we define, the constant
5.15we 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.17Consider the instance when
5.18where Re[·] denotes taking the real part. Equation (5.18) implies that An integral relation, only true for and is
5.19where *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.20For 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.21For *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.22On the other hand, for *x* < 0, we use the contour in figure 11*c* to enclose the pole This results in
5.23Combining 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}*.

## 6. Illustrative applications

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

*ω*(

_{j}*k*). The components

*T*are then used to analytically derive the Green function component

_{j}*G*, both in two dimensions and in one dimension, that will significantly reduce the complexity of calculating the BOLD response to any neural drive. In this section, we illustrate how our theoretical derivations can be used to calculate the components of the haemodynamic response to a specific neural stimulus; in this case, we use a spatio-temporal Gaussian neural stimulus. We also show an example neural stimulus, a spatio-temporally varying sinusoidal stimulus, that can potentially enhance a BOLD response component of interest. For illustration purposes, we will restrict our analysis to one-dimensional haemodynamics.

_{j}### 6.1. Haemodynamic response to a spatio-temporal Gaussian neural stimulus

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

*Γ*to predict the BOLD response of one subject and compare with their data. We will then study the properties of the resulting decomposed BOLD components.

We model the stimulus mathematically as
6.1where *σ _{x}* and

*σ*

_{t}are the spatial and temporal widths, respectively, and

*t*

_{0}is a temporal parameter to match the neural stimulus' centre. In this case, we use and to have fixed full widths at half maximum of 2 mm and 2 s, and fix

*t*

_{0}= 4 s, to emulate the neural response elicited by the experiment of [13]. We then convolve equation (6.1) with the one-dimensional Green functions (equations (5.12) and (5.24)), using

*v*= 2.01 mm s

_{β}^{−1}and

*Γ*= 0.72 s

^{−1}, to get the responses for one subject as shown in figure 12. Note that the chosen values of the parameters

*v*and

_{β}*Γ*are taken from the experimental estimates of [13].

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 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.

### 6.2. Enhancement of blood oxygen-level dependent response components

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.2where *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}*(

*x*,

*t*) or substitute equation (6.3) into equation (3.23), resulting in 6.4We further simplify equation (6.4) by substituting equation (3.2) and exploiting the fact that

*a*(

_{j}*k*) and

*ω*(

_{j}*k*) are even functions, to get 6.5where

*a*(

_{j}*k*) is defined in equation (3.17). Equation (6.5) can be simplified to 6.6and 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}*(

*x*,

*t*). First, we can vary the spatial frequency

*k*

_{s}of the neural stimulus such that it is tuned to the resonance of

*a*(

_{j}*k*). That is, we use where is the analytical resonance frequency we derived in equations (4.8)–(4.10). However, for pairs of

*v*and

_{β}*Γ*where becomes imaginary and this technique will not be applicable. Second, we can manipulate the temporal frequency

*ω*

_{s}to minimize the denominator of equation (6.6), such that However, because

*ω*

_{s}is real-valued, it can only have as its limiting value.

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.8which 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.

## 7. Summary and conclusions

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

*Γ*by exploring the wave components of the BOLD response and disregarding the non-wave components.

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

*Y*. Exploration of this application and possible nonlinear BOLD dynamics is left for future work.

_{j}## Authors' contributions

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.

## Competing interests

The authors declare no competing interests.

## Funding

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).

## Acknowledgements

The authors thank Thomas Lacy for useful discussions.

- Received March 31, 2016.
- Accepted April 12, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.