## Abstract

It is shown that recently discovered haemodynamic waves can form shock-like fronts when driven by stimuli that excite the cortex in a patch that moves faster than the haemodynamic wave velocity. If stimuli are chosen in order to induce shock-like behaviour, the resulting blood oxygen level-dependent (BOLD) response is enhanced, thereby improving the signal to noise ratio of measurements made with functional magnetic resonance imaging. A spatio-temporal haemodynamic model is extended to calculate the BOLD response and determine the main properties of waves induced by moving stimuli. From this, the optimal conditions for stimulating shock-like responses are determined, and ways of inducing these responses in experiments are demonstrated in a pilot study.

## 1. Introduction

Functional magnetic resonance imaging (fMRI) is commonly used to examine changes in neural activity by detecting changes of the blood oxygen level-dependent (BOLD) signal in brain tissue due to stimuli. In order to better understand the mechanisms that are important in the generation of the BOLD response, quantitative physiological-based modelling is thus required. Early models were restricted to modelling the temporal structure of the response [1] while assuming that the spatial component was separable from this. Recent advances in modelling along with newer experiments have shown that the haemodynamic response is coupled in both space and time [2–5]. This is due to blood flows between adjacent regions in the brain resulting in spatial couplings; in particular, a new theory indicated that haemodynamic waves could be induced by appropriate stimuli, as was confirmed by experiments [3].

A recent model [2,6] that includes these effects obtained the BOLD response to stimuli by using mean field (small-scale effects are averaged over) and poroelastic (modelling the cortex as a porous medium perfused by blood) approximations to derive a spatio-temporal model for the BOLD response. By using appropriate equations for arterial inflows and deoxyhaemoglobin (dHb) dynamics, as well as applying conservation of mass, momentum and oxygen, a system of differential equations was derived that can be solved to predict the response. The linear response function of the system implies that haemodynamic waves can be induced within the cortex, and that these waves propagate several millimetres with a velocity of the order of 1 mm s^{–1} before becoming too small to be detected due to damping. These predictions were subsequently validated by experiment [3].

The existence of haemodynamic waves immediately implies that moving the neural stimulus will tend to induce ‘piling up’ of wavefronts as happens with bow waves of ships in water and sonic booms produced by aircraft. This effect is expected to be strongest when the velocity of the drive across the cortex closely matches or exceeds the velocity of the waves generated by the stimulus, as the front of the wave will always be at the same location as the neural stimulus, and any waves generated by the stimulus can directly reinforce the wave. This type of stimulus should thus produce a significantly stronger BOLD response than other stimuli and could be intense enough to demonstrate nonlinear effects predicted by the theory.

This paper uses haemodynamic theory to predict induced BOLD responses without making linear approximations, because some papers have reported that the nonlinearities are significant [7], while others have disagreed and reported that the response is approximately linear [3]. This means that the threshold where nonlinear effects become significant should be determined to understand when the assumptions used in linearizing are inaccurate. If these nonlinearities are sufficiently strong, a nonlinear steepening of the front of the travelling wave is also possible, resulting in a shock-like haemodynamic response (has a very steep wavefront). The peak of the wave should also saturate nonlinearly at some level, because there is a limit to how much the blood vessels in the brain can expand. By analogy with shocks that occur with waves in other mediums, forming a shock would require the stimulus to be reasonably strong and very localized, while being moved at or above the haemodynamic velocity. These constraints mean that if a shock-like response can be induced via nonlinear interactions, the precise shape of the shock and the conditions under which it can form will depend on the physical parameters of the brain to which the neural stimulus is being applied, thereby allowing for the parameters to be measured non-invasively. Use of modelling to optimize the applied stimuli has been done before, but only for the temporal aspects of the BOLD response [8].

Moving neural drives and stimuli are central to experimental inference of retinotopic mappings of the brain [9], and determining population receptive field (pRF) sizes in the visual cortex [10], which involve mapping of the visual field onto the visual cortex. Thus, it is possible that haemodynamic waves are already being exploited when moving stimuli are used to increase the induced BOLD response and improve measurements of activity in these experiments. The moving stimuli generate significantly stronger waves than can be obtained for stationary stimuli, and suggest the quality of the data in retinotopic mappings and pRF modelling has been improved by inducing spatio-temporal resonances. This is because stimuli that are moved at the same velocity as the waves they produce will have their response especially reinforced by these waves.

In this paper, we show how the nonlinear equations of the model can be solved in order to obtain the predicted BOLD response for a given neural drive. Then we predict the responses for certain physiologically relevant ranges of physical parameters, and determine under what conditions the nonlinear effects will produce the strongest responses to moving neural drives. Finally, we use the results to improve experimental designs, and compare them with current experiments.

## 2. Methods and theory

In this section, the physical principles used to derive the haemodynamic model linking visual stimuli with the measured BOLD response are stated, along with the equations derived from them, and the method used to perform calculations using the model is outlined.

### 2.1. Haemodynamic model

Here we outline the model used, whose detailed derivation and discussion have been presented elsewhere [2,3], explain briefly how each equation is derived from physiology and link them together as shown in figure 1.

Neural activity can drive haemodynamic activity. Neural activity, specifically the increase of neural firing and synaptic activity, can induce changes in the local arterial inflows of blood, due to neurovascular coupling between the neurons and the blood vessels via astrocytes and neurotransmitters. The signal that causes vasodilation in response to neural activity can be modelled as a neural drive. However, because the precise dynamics of the link between the neural activity and the neural drive are not fully understood yet, we start with a given profile for the neural drive itself, and then determine the BOLD response *Y*(**R**, *t*) that can be measured in an fMRI scanner if a stimulus that produces that drive is applied. The drive results from the effect of dynamics in such a way that it can be approximated as the neural activity, with spatial spreading and a low pass filter applied. This means we can define an appropriate neural drive that can be experimentally induced.

We consider a block of cortical tissue as seen in figure 2. The main blood inflows occur at a depth *z*_{0} ≈ 1.3 mm below the surface, as this is the layer with the highest blood mass density and the local demand for oxygen is assumed to be proportional to the blood mass density in the layer. Note that although the inflow arrives from descending arteries, the modulation to the mean field of blood mass occurs where the main source of vasoactive agents occur. The cortical tissue is approximated as a poroelastic medium perfused by blood, as shown in figure 3. This means that the local concentration of blood within the cortex can be related to the local pressure. In addition, small-scale effects (smaller than approx. 1 mm in linear scale) are averaged over in this model, which thus makes a mean field approximation.

The arterial blood inflow rate [11] *F*(**R**, *t*) is modelled as having a damped harmonic response driven by the neural drive , with [2,12,13]
2.1where *κ*, *γ* and *F*_{0} are the blood flow signal decay rate, the flow-dependent elimination constant and the arterial blood inflow with no applied stimulus, respectively.

The total mass density of the brain is the sum of contributions due to the brain tissue and blood. The mass density of blood in the cortical tissue is calculated by relating it to the blood pressure, using conservation of the mass and momentum of the blood within the tissue. A key step is relating the pressure within the blood vessels *P*(**R**, *t*) to , which is done via the constitutive equation
2.2where *β* is the elasticity exponent of the blood vessels and *c*_{2} is an empirically derived constant of proportionality. This gives the pressure in terms of blood mass density, assuming that the tissue is a continuum. A value of *β* = 1 corresponds to the vessels being perfectly elastic, while the value of *β* ≈ 3.23 used here (from the inverse of Grubb's exponent, which relates cerebral blood flow to the cerebral blood volume) implies that they are hyperelastic (i.e. more resistant to stretching than perfectly elastic vessels). The value of *β* depends on the type of vessel being considered (arterioles are more elastic than veins, and so have a lower *β*), but as this is a mean field model, we take an average value for all vessels.

Conservation of blood mass within the brain implies [6]
2.3which relates the rate of change of the blood mass density within the cortical tissue to inflows () and outflows () of blood (*ρ _{f}* is the blood mass density with no tissue present), and where

**v**(

**R**,

*t*) is the average velocity of the combined tissue and blood. The function

*g*(

*z*) has been set to be a Gaussian centred at

*z*

_{0}= 1.3 mm, to reflect the higher level of inflows in the layers near that depth, with full width at half maximum (FWHM) of approximately 2.4 mm (the inflows near

*z*= 0 mm and

*z*= 3 mm should still be significant, which is why the FWHM is comparable with the cortical thickness). The dependence of at equilibrium with respect to

*z*was determined by measuring the total vascular volume [14] and then finding a profile that matches the measured values.

Momentum is also conserved for the blood. By equating the forces acting on the blood with the rate of change of momentum density, one obtains
2.4where *c*_{1}/*ρ _{f}* is the constant of proportionality between blood acceleration and pressure gradients in the medium and

*D*quantifies damping due to blood viscosity. In the case of steady blood flows, the left side is zero and this equation reduces to Darcy's law from poromechanics [15].

From equations (2.3) and (2.4), one obtains the following nonlinear equation for [3]:
2.5which relates to *F*(**R**, *t*), because *P*(**R**, *t*) is a function of via equation (2.2).

The concentration of dHb, *Q*(**R**, *t*) can be determined by relating its time derivative to its flux into nearby regions, its rate of increase from the consumption of oxyhaemoglobin (oHb) by neurons, and its rate of decrease due to outflows of blood. This gives an equation that relates *Q*(**R**, *t*) to and the boundary conditions of the system, with the rate of flow out of the cortical tissue determining how fast the dHb is cleared from the brain. So, haemoglobin conservation yields
2.6with *ψ* being the ratio of haemoglobin concentration to blood density in mol kg^{–1} in SI units and *η* being the fractional oxygen consumption rate per unit time in s^{–1} in SI units. Recent works [16,17] indicated that vein diameters do not change significantly during stimulations of less than approximately 30 s, compared with changes in arteriole diameters. Hence, we assume that the changes in blood mass density (as a contribution to total density) are primarily confined to the arterioles, while the changes in dHb content are confined to the veins (as the main transfer of oxygen from the blood to the neurons occurs in the capillaries). This means that the divergence term, as needed in a continuity/conservation equation, corresponding to changes in dHb due to blood flow to adjacent areas of the brain, can be neglected, because that blood flow will contain little dHb. This is a change in the approximations relative to our previous papers [3,6] in which the divergence was treated as being very small, but potentially non-zero. Comparing the results simulated with and without the divergence term from previous papers shows that there is very little difference in the responses induced, so it is omitted henceforth.

Finally, the BOLD signal *Y*(**R**, *t*) is determined from and *Q*(**R**, *t*) using a semiempirical relation derived from the properties of the fMRI scanner such as the magnetic field strength and the method of signal acquisition [18]. This gives
where *V*_{0} is the resting blood volume fraction, *Q*_{0} is the resting dHb content per unit volume and *k*_{1}, *k*_{2} and *k*_{3} are parameters that depend on the fMRI scanner being used and the experimental protocol [18,19], as this changes how sensitive it is to intravascular signals compared with extravascular signals.

Using equations (2.1)–(2.7), it is possible to predict the measured BOLD response *Y*(**R**, *t*) to a given neural drive by calculating the intermediate quantities *F*(**R**, *t*), , *Q*(**R**, *t*) and then substituting them into equation (2.7). The actual methods of calculating the nonlinear response are detailed in appendix A, including the specific finite differencing method used to solve the differential equations. In what follows, all of the values of constants involved are either obtained from best known physiological estimates or indirectly from these by solving our equations for a steady state where no neural drive is applied, and all time derivatives are thus set to zero [6]. These values appear in table 1.

### 2.2. Moving stimuli

In order to induce haemodynamic shock-like waves, we consider stimuli that cause a drive that moves along the surface of the cortex. For the purposes of inducing the strongest shock-like response, we consider line stimuli moving across the cortex, as shown in figure 4, because point stimuli will induce BOLD responses that tend to spread out in all directions [6], whereas for line stimuli, spreading can only be perpendicular to the line. Without the reinforcement provided by adjacent points with stimuli moving in the same direction, there would be a smaller induced response and smaller nonlinearities in the case of a point source.

By comparing the situation with waves in other media, we predict that the shape of the response should depend significantly on the ratio of the speed of the neural drive across the cortex to the haemodynamic wave velocity. The strongest response should occur when the drive is moving at approximately the same velocity as the waves induced by the drive, thereby inducing a resonant response. With planes travelling through air, this corresponds to a plane moving at Mach 1, so we introduce a haemodynamic Mach number *M*_{H} to describe the velocity of a moving drive, defined by
2.8and determine how the response varies with *M*_{H}.

In order to induce moving line drives in cortex, we make use of the well-known retinotopic mapping from the visual field to cortex [25]. These maps show that a circle, placed at the centre of the visual field, projects to an approximately straight line in primary visual cortex. These relationships have been quantified by various studies [26], leading to the analytic expression
2.9where *w* is the location in V1, *k* is a scaling factor, *E* is the retinal eccentricity (in degrees), *P* is the angle between the point in visual space and the horizontal (in degrees), *f _{a}* and

*f*are shear functions from the double-sech model of Schira

_{b}*et al.*[26] in degrees

^{–1}, and

*a*and

*b*are parameters that define where the foveal singularities are (at 1.05° and 90°, respectively). This means that if expanding or contracting toroidal sectors (arc segments) are presented, as in figure 5, they will each induce a moving line drive. We can use equation (2.9) to estimate the speed of expansion (or contraction) needed to generate a neural drive travelling at a constant speed

*v*. It is thus possible to determine the visual stimulus required to generate arbitrary neural drives, but we only examine point and line drives in this paper.

_{n}The typical drive used in the simulations is a spatial Gaussian, which is applied for an interval of length *τ*_{s}. Explicitly, it is described by
2.10where *v*_{N} is the drive's velocity across the cortex, *σ* is the spread of the Gaussian (the FWHM is approximately 2.35 times this and was set to 1 mm) and **R** = (*x*, *y*, *z*) is the location of the point being considered in three-dimensional space. In the case that a one-dimensional stimulus is being modelled, the stimulus is independent of *y*, so the drive used for bar stimuli is
2.11

Because a point in the visual field induces a drive that is approximately Gaussian with a FWHM of 1 mm, this is the smallest spreading that the response can have. The amplitude of the drive is chosen so that the largest change in *F*(**R**, *t*) is a doubling of *F*(**R**, *t*), which is approximately the maximum normally observed in experiments [27]. This type of stimulus is similar to the ones used in retinotopic mappings [9].

### 2.3. General magnetic resonance imaging procedures

For the imaging data, we performed a pilot study with one healthy subject (male aged 25) in order to illustrate the effects of our findings. By restricting the fMRI scans to the occipital pole, we achieved a resolution of 1.5 × 1.5 × 1.5 mm^{3} and 2 s TR echoplanar images (EPI). These EPI data were then coregistered to a high-resolution T1-weighted anatomical scan, acquired at 0.75 × 0.75 × 0.75 mm^{3}. Furthermore, these mappings were restricted to the grey matter by segmenting the anatomical data into grey and white matter. Finally, functional retinotopic scans were used to map out the expected cortical positions in the visual cortex of the subject. Data were acquired on a Philips 3 T Achieva Series MRI machine equipped with Quasar Dual gradient system and a 32-channel head coil.

### 2.4. Visual paradigm

The stimuli consisted of expanding wedges, as seen in figure 5, that were made to induce neural drives that travelled at various speeds on cortex. We aimed to produce neural drives that were moving bars of width 3 mm and height 10 mm that travelled 30 mm. To this, we produced a series of 50 bars where the centre of each bar was placed 0.6 mm from the centre of the previous bar. By using equation (2.9), each bar corresponds to a wedge in visual space, where the stimulus started: at 0.3° eccentricity, spanning 15° polar angle, and ended: at 5.5° eccentricity, spanning 7° of polar angle. Presenting each of these wedges sequentially at 1200 ms corresponds to neural drives travelling at 0.5 mm s^{–1}. In order to increase activity, we added a checkerboard pattern which alternates, by contrast, every 200 ms, and we also included static concentric rings and a simple fixation task which requires the subject to press a button when the 4 × 4 pixel square in the centre changes to red. The stimuli are synchronized with the scan acquisitions to start at the first volume, with a 20 s rest period and are repeated three times, thus lasting 240 s.

This procedure is repeated for neural drive speeds of 1, 2, 3, 5 and 10 mm s^{−1} with 5, 7, 8, 9 and 10 repetitions lasting 250 s, 245 s, 240 s, 234 s and 230 s, respectively. The stimuli are coded in Matlab with Psychtoolbox (http://psychtoolbox.org/) and presented on a screen which had a range of 11° of eccentricity.

### 2.5. Anatomical data and preprocessing

A high-resolution anatomical T1-weighted image is acquired over the entire brain using a MPRAGE imaging sequence [28]. We acquire 250 0.75 mm slices, with a 340 × 340 matrix and a 256 mm field of view, resulting in an isotropic resolution of 0.75 mm in each direction. This image is automatically segmented according to pixel intensity, guided by anatomical templates using Freesurfer (http://freesurfer.net/). Cortical surfaces, for the left and right hemispheres, are reconstructed from the segmentation. The first is at the interface between the grey and white matter denoted by ‘white’. The second represents the pial surface and is generated by expanding the white surface to closely match the interface between the grey matter and cerebrospinal fluid. Topological defects in these two surfaces are corrected, then smoothed to closely follow the two interfaces. One additional surface is constructed by expanding the white surfaces by 50% of the cortical thickness along surface normals, denoted ‘midlayer’.

Two flattened maps are constructed around the occipital pole of the right and left hemispheres by taking a patch of 80 mm from the fovea of the left and right hemisphere, and projecting the vertices onto a circle using Tutte's embedding theorem [29], an improved version of the procedure used in Aquino *et al.* [3].

### 2.6. Functional data and preprocessing

In the fMRI acquisition, we use an accelerated EPI sequence, using a SENSE (sensitivity encoding) factor of 2.3 [30] on one subject. Functional data are acquired in 32 1.5 mm slices, with a 128 × 128 matrix and a 192 mm field of view. Our functional scans are acquired at a time of repetition of 2 s, and an echo time of 28 ms.

Functional data are motion corrected and slice scan-time corrected using mrTools (http://www.cns.nyu.edu/heegerlab/wiki/doku.php?id=mrtools:top). The data are then averaged according to the number of repetitions allowed by the visual paradigm. This means that the results for 0.5 mm s^{–1} are the average of three repetitions, while the results for 1 mm s^{–1} resulted from the average from five repetitions due to each repetition taking less time, and higher drive velocities allowed for even more repetitions. We reconcile the discrepancy between the averaging conditions in that the neural signal will be stronger for lower velocities thus higher signal-to-noise ratio requiring fewer averages and thus the inverse for the higher velocities. In addition, after this averaging we low pass the time series using a Butterworth filter of order 3 and low-pass frequency of 0.1 Hz.

The fMRI data are analysed in three different spaces. (i) The original planar space of the data acquisition. (ii) The cortical surface—defined by the middle grey layer (see above) into which the data were interpolated, resulting in spatio-temporal data on the left and right hemisphere cortical surfaces. We note that this required alignment of the EPI data onto the high-resolution T1 scans where the surfaces are generated and this ensures only one interpolation step onto the cortical surface directly from the EPI data. (iii) The last of these were on the averaged coordinate space along the direction of the neural drive. This is achieved by taking the visual stimulus and finding the polar angle and eccentricity of each stimulus wedge. As the induced BOLD response of the stimulus projects past the retinotopic location of the stimulus, we extend the mapping to find corresponding retinotopic locations if the stimulus extends to 15° of eccentricity as a way to find an averaging region that respects the organization of visual cortex. The corresponding locations on cortex are determined by making use of a retinotopic atlas [31] which estimates eccentricity and polar angle on the cortical surface from folding patterns in striate cortex. Each stimulus wedge is then mapped onto a small strip resulting in a series of parallel lines on cortex.

The time series from each of these lines are averaged, resulting in a mapping from the surface to one dimension in space. The shortest geodesic distance between each line and the first line (defined to be the one closest to the fovea) defines the distance coordinate *x*. No further preprocessing of the data in the spatial dimension is performed. We stress that the mapping between the cortical space and the flattened space is a one to one mapping purely done for visualization purposes.

## 3. Results

In this section, we predict and compare BOLD responses for stationary and moving neural drives, in both the linear and nonlinear cases of the model, to determine which produce the responses that provide the most useful information. In addition, we determine the situations in which the strongest shock-like responses are induced, what sort of drive needs to be applied in order to produce them, and their properties.

### 3.1. Types of drives

The first kind of drive we consider in figure 6 is a moving bar drive (as was depicted in figure 4), so we can examine the BOLD response induced by such a drive when it is moved. Using this kind of drive gives the set of responses shown in figure 6, with each row showing a different value for *M*_{H} going from 0 to 1.5. When *M*_{H} = 0.5, i.e. on the second row of figure 6, the BOLD response starts off at baseline, then increases in amplitude while moving to the right, before plateauing while still moving. After about 9 s, the profile of the response stays constant, but moves to the right, with a negative wake trailing behind it. This overall structure is similar for all moving stimuli; however, for a stationary stimulus (*M*_{H} = 0) we only see the spreading of the response due to haemodynamic waves. The equilibrium response that the *M*_{H} = 0 drive reaches after about 9 s is smaller than the moving responses for *M*_{H} = 0.5 and 1. As seen in all these responses, the BOLD response is approximately uniform in the zone −5 mm < *y* < 5 mm.

The constant profile with respect to *y* suggests that if the length of the bar is sufficiently large, then the *y* direction is irrelevant for parts of the response near the centre of the bar. In this case, the response can be modelled using the one-dimensional drive in equation (2.11) instead of a two-dimensional one. Because this gives one less dimension for the haemodynamic waves to propagate in, it leads to less spreading of the response and thus stronger BOLD responses and nonlinear effects. We next examine the space–time plot for the one-dimensional responses with the same parameters as in figure 6 to determine the dynamics. Figure 7 shows that the induced response starts off increasing, before moving significantly to the right. The induced response reaches a moving steady state after approximately 10 s, and the response dies away after it is removed at 15 s. For *M*_{H} = 0 and 0.5, figure 7*a,b*, the only visible negative BOLD response occurs after the drive has been removed, and it takes about 3 s after the removal for the response to drop below 0. This delay is typical for BOLD responses in general. However, it can be seen that in the cases of *M*_{H} = 1 and 1.5, in figure 7*c,d*, a negative response trails the main wavefront, which is moving to the right. This is due to the negative induced response being left behind by the moving drive, thereby allowing it to be seen seperately. In the *M*_{H} = 0 and 0.5 plots, in figure 7*a,b*, the negative part of the induced BOLD response significantly overlaps the positive part, which masks it. This also explains why the response for figure 7*a* peaks at around 6 s because beyond that time the positive response induced by the drive starts overlapping with the negative response induced by the drive at earlier times, giving a lower equilibrium response until the drive is removed at *t* = 15 s.

To see that the response does reach a moving steady state after 10 s, we can also shift into the reference frame of the drive, by replacing the *x* coordinate by . The resulting responses are shown in figure 8. In this frame, it is seen that the response reaches a steady state after about 10 s, and remains in this state until the drive is removed. It is also important to note that the response is not able to extend very far into positive *w* values when *M*_{H} > 1 in figure 7*d*, because the drive is moving as fast as or faster than the waves; the response at *w* > 0 is due to the drive having a FWHM of 1 mm.

### 3.2. Effects of changing drive velocity

To further examine the effects of changing *M*_{H} beyond those already seen, we note that figure 8 shows that the response to a moving drive eventually reaches a steady state in the frame of reference of the drive. Thus, we can examine how the cross sections of the steady-state moving responses change with *M*_{H}. If certain values of *M*_{H} give a particularly strong response, then moving a drive at that velocity gives a method for inducing a particularly strong signal that can be more easily detected against noise in experiments.

Figure 9 shows the profile of the BOLD response to a one-dimensional line stimulus moving at a constant velocity after it has reached a steady state in its comoving frame of reference for different values of *M*_{H}. As *M*_{H} increases, the profile initially narrows due to the movement of the drive not allowing the spatial spreading to completely stabilize. However, as *M*_{H} increases further, the response becomes more stretched out along the *x* direction, and becomes more and more asymmetrical, with an increasingly large tail for negative *w* compared with the positive tail. This stretching is due to the drive moving faster through the tissue, so any particular point is stimulated for a shorter time and thus has a smaller induced response.

The strongest BOLD response seen in figure 9 is for *M*_{H} ≈ 0.75. This is due to how the two separate mechanisms that limit the height of these profiles change with *M*_{H}. The first mechanism involves the fact that a stationary stimulus will produce a negative response shortly after it is removed (as seen for *M*_{H} = 0 in figures 7 and 8). This means that if a drive is continuously applied to one location, the negative response induced at earlier times will overlap with the positive response induced later, thereby making the equilibrium response lower overall (seen again in the same figures, as the response peaks after approx. 5 s and then decreases even though the drive is still being applied). The second mechanism relates to the fact that the faster the stimulus is moved, the more spread out the response is in the direction of motion as mentioned before, resulting in smaller peaks. The combination of these effects means that the greatest steady-state response should occur for a value of *M*_{H} between 0 and 1 because if *M*_{H} is greater than 1 then the wavefront will always lag behind the location of the drive, decreasing its size. For the parameters used here the peak response is at *M*_{H} ≈ 0.75, as seen in figure 9, and is in accord with the above arguments.

### 3.3. Strength of nonlinear effects

In order to induce shock-like responses, the drive should be as strongly localized in the direction of travel as possible, while extending as far as possible in the perpendicular directions. This is to ensure that the sharpest wavefront can be obtained, due to the localization, while the extension in the other dimension causes it to approximate a line stimulus, and hence makes the induced response stronger.

Changes to tissue and blood parameters can also change the strength of nonlinearities that are present. Low viscous damping (low *D*) allows for larger changes in blood mass density than when there is high damping, because high damping causes the haemodynamic waves to decrease in amplitude faster. The drive should be moved at slightly above the equilibrium haemodynamic velocity *v _{β}* when attempting to induce a haemodynamic shock, because the nonlinear pressure term increases the haemodynamic velocity in areas of increased blood mass density . This means that the front of the haemodynamic wave moves slightly slower than its crest, thereby leading to nonlinear steepening.

The size of these nonlinear steepening and damping effects can be directly tested by comparing the responses predicted for a drive with and without the nonlinearities present. The nonlinear effects can be effectively removed relative to linear effects by decreasing the amplitude of the drive by a factor of 1000. Because the nonlinear effects are roughly proportional to the drive amplitude to the power of *β*, this reduces them by a factor of roughly a billion. Multiplying the final predicted response by 1000 again restores the linear effects to their normal strength, but the nonlinear effects are still reduced by a factor of a million, giving a numerical linearization of the equations (linearizing the equations explicitly can also be done, and was done by Aquino *et al.* [3]). Figure 10 shows that when the drive is moved slightly above the haemodynamic velocity, the main effect of nonlinearity in the model is to cause a saturation of the response at a lower level than the linearized version. This reduction in the response is significant (by approx. 30% of the linear response), and opposes steepening of wavefronts, and hence shock formation. The difference in the responses shows the linear response is more positive than the nonlinear response everywhere except the very front of the wave, where the nonlinear response is slightly more positive. The ratio of the responses, shown in figure 10*b*, indicates that the magnitude of the nonlinear response is slightly higher near the very front of the wave. However, it then steadily decreases moving towards the crest, meaning that the response is slightly steeper for that case, and that the size of the wake behind the wave is always at least 10% larger for the nonlinear response. For shocks to exist in the model, the effects of the nonlinear steepening would have to be stronger than the saturating effects near the crest, to give a steeper wavefront than the linear model, but figure 10 indicates this is not the case.

To understand more why this flattening of the nonlinear response is occurring, we examine the profiles of and the pressure, because the relationship between and pressure is the main nonlinearity in the system of equations. These profiles are shown in figure 11, and shows that while the values follow approximately the same curves as the BOLD responses do in figure 10, the pressure in both cases is almost identical. This suggests that the difference between the linear and nonlinear terms can be estimated from the percentage change in pressure. In the linear case, we have 3.1and hence 3.2

In the nonlinear case, we have equation (2.2) or equivalently, 3.3so the change for the nonlinear is 3.4

If we have a stimulus that gives an 87% increase in pressure, as in figure 11, then the above equations imply that the increases in will be 27% and 21%, respectively, as also seen in figure 11. This means that there is a 25% larger increase for the linear case, as is the case for the stimulus in figure 10, and hence we can estimate how large the discrepancy between the two cases will be from the magnitude of the drive and from the value of *β*.

We can verify that the steepening will not be able to produce shocks even under extreme stimuli by reducing the damping from 200 to 50 kg m^{–3} s^{–1} to increase the response, and also replacing the drive with
3.5where *H*(*t*) is the Heaviside step function (equal to 1 if *t* ≥ 0, 0 otherwise). This drive is a step moving to the right at velocity *v*, and with a discontinuity at *w* = 0 which should give the steepest possible wavefront within the model. Using this drive gives the response shown in figure 12, which demonstrates that even with a discontinuity at the front edge of the moving drive, and very low damping, the response is unable to steepen to the point at which a shock front is formed. The failure of these highly favourable conditions to result in a shock indicates that shocks almost certainly cannot be produced, subject to our model equations. However, as the steepest wavefront has a gradient of approximately 1% change in BOLD signal per millimetre, the wavefronts of these responses would be very easy to detect, so attempting to reproduce this kind of drive would still be useful.

### 3.4. Parameter dependences

It is also important to examine the effects of changing parameters within the model, especially when these parameters vary significantly between subjects. The parameter most important for determining the shape of the response induced that varies significantly between subjects is *D*, the effective blood viscosity, which is related to the viscosity of blood *μ* and the permeability of the cortical tissue *λ* by . In addition, *D* can be changed by drugs that change *μ*, such as pentoxifylline and oxypentifylline, which could potentially allow experiments to be done on the same subject for different values of *D*.

Figure 13 shows how changing the value of *D* changes the profile of a moving wavefront in the reference frame of the moving drive, with *M*_{H} = 1. The shape of the main peak in the response is nearly constant, though slightly higher and further forward when *D* is lower. The main effect of changing *D* is that the size of the negative response behind the drive changes significantly, with the size being largest for low *D*, and a much smaller wake for higher *D*. This indicates that in most situations the wake will not be visible, and would require low blood viscosity to be observed.

### 3.5. Experimental data

We can compare the predictions from the previous simulations with experimental data obtained by using the visual stimuli shown in figure 5 to induce moving bar drives in a subject. The value of *M*_{H} (*ν _{β}* is approx. 2 mm s

^{–1}for this subject, from a previous experiment involving this subject [3], so

*M*

_{H}= 1 corresponds to

*ν*

_{N}= 2 mm s

^{−1}) is changed by changing the rate of expansion of the stimulus ‘wedges’ in the visual field, which changes the rate at which the bar drive moved across the cortex. The measured responses are shown in figure 14, shown in the same form as in figure 7. This shows that the induced response is reasonably strong for low drive velocities (first panel with

*M*

_{H}= 0.25), increases slightly for

*M*

_{H}= 0.5 and then steadily decreases as the value of

*M*

_{H}continues to increase, which agrees with the predictions from figure 9 up to

*M*

_{H}= 2. When

*M*

_{H}is increased beyond 2 in the model, the overall shape of the response remains constant but becomes more stretched out in the

*x*direction, with a lower peak response, which matches what is seen here in the last three panels, as they all have a very similar overall shape, with the main difference being a lower BOLD response for faster drives. This demonstrates that there is little point in moving a drive at speeds greatly above the haemodynamic velocity, as in those cases the BOLD response induced will be very small and difficult to measure, compared with the case of lower

*M*

_{H}.

There is also a significant negative BOLD response present in all of the panels, especially after the drive has passed by a section of the cortex. This matches the predictions made in figure 7, with the negative response peaking at a magnitude of about 30% as high as the main peak of the positive response. This is consistent across all of the values of *M*_{H}, and indicates that these kinds of stimuli are also good for inducing strong negative responses.

## 4. Summary and discussion

A spatio-temporal model for the induced BOLD response to moving stimuli was solved numerically including nonlinear effects, particularly looking for responses that correspond to phenomena similar to shock waves, and the resulting responses were analysed. The main results are:

(i) Moving a neural drive shaped like a bar in V1 perpendicular to its length, as in figure 4, results in the BOLD responses shown in figure 6. These show a large positive wave moving slightly behind the bar with a negative wake behind the wave. Positive responses spread outwards at the haemodynamic velocity

*ν*to the sides beyond where the drive was applied. For stimuli that are not bar shaped, these waves still exist, so determining their properties is important for improving fMRI studies in all regions of the brain._{β}(ii) The BOLD response for a moving bar drive is translationally invariant along the bar's length except near its ends. In this case, only one slice of the drive and response has to be modelled, because it reduces to a one-dimensional problem, with the behaviour being constant along the bar's length.

(iii) For a given set of parameters, a family of response profiles has been produced by moving a neural drive across V1 at different drive velocities,

*ν*_{N}, as shown in figure 9. This demonstrates that for prolonged stimuli, the strongest response is not for a stationary stimulus, but for one moving across the cortex at (for the parameters of figure 9, it was strongest for ). This dependence gives a way of determining*ν*and other related properties non-invasively (this has also been done by measuring the angle at which the wake spreads outwards, as in Aquino_{β}*et al.*[3]), by applying drives for different*ν*_{N}to reproduce the response profiles shown in figure 9. From this one can then determine*ν*. This dependence can also aid experimental design by suggesting improved visual stimuli optimized to individual subjects, with_{β}*ν*_{N}so that responses are either minimized or maximized, depending on what is preferred in a given experiment.(iv) Nonlinear effects result in lower amplitudes for the induced BOLD response than those predicted from linearized equations. This also applies to the underlying quantities (such as the local blood mass density and the dHb concentration), indicating that the nonlinearities tend to limit the amplitude of the BOLD response and will not produce shocks. This is because pressure terms in equation (2.5) are larger when not linearized, and increased local pressure tends to force blood into adjacent areas or out of the brain through veins, decreasing the local blood mass density and BOLD response. This means that true haemodynamic shocks probably do not exist; indeed, if nonlinearities had instead tended to strengthen the response they would have probably been seen experimentally by now.

(v) Comparing the model predictions in figure 7 and the experimental results for bar drives in figure 14 shows that the model reproduces the main features of the BOLD response morphology and amplitude for several different

*ν*_{N}.(vi) Our numerical calculation of the response includes the direction perpendicular to the cortical sheet, unlike the

*z*averaging that has been done previously and demonstrates that earlier assumptions that the*z*-direction can be dealt with by assuming separation of variables is not completely true. This is because the waves in the blood mass density and other quantities being calculated cannot propagate instantaneously through the 3 mm of cortex in that direction, but instead travel at*ν*≈ 1 mm s_{β}^{−1}. This means that treating the*z*-direction as separable will result in inaccuracies, because the waves cannot propagate instantly through the layers.

Future work with the existing model will require more moving-drive data for a larger number of subjects and with a range of values of *ν*_{N} close to *ν _{β}* for the subject (values of do not yield much useful information because the BOLD signal amplitude becomes too low). With these data available, the procedure of determining

*ν*in the subjects given the moving stimulus data can be tested properly, in addition to the predictions resulting from figure 9 regarding the values of

_{β}*ν*

_{N}that lead to a maximal BOLD response.

In the future, there are several ways in which the model can be extended and improved. Hysteresis effects for the pressure–volume relationship could be included to more accurately reflect how the cortical blood vessels respond to large time-dependent variations in the blood mass density. The modelling of arterial inflows could be made more detailed instead of just assuming that the inflows can be approximated as following damped harmonic oscillations, and by modelling the behaviour of astrocytes in more detail. Finally, if smaller scale effects (less than or equal to 1 mm) are to be modelled, it will become increasingly important for the averaging that is implicit in the mean field approximation to be modified to take into account local variations in the types of blood vessels present in the tissue. Such variations in the proportions of blood mass within arterioles, capillaries and veins can change the BOLD response measured significantly, due to the different levels of oHb and dHb in each kind of vessel. This is especially important for veins because a large vein can show up as having an inverted BOLD response compared with the surrounding tissue.

## Ethics

The study protocols were approved by ethics boards of the University of New South Wales and Neuroscience Research Australia.

## Authors' contributions

The idea was conceived by P.A.R. The modelling and analysis were carried out by T.C.L with supervision from K.M.A and P.A.R. The experiment was designed, carried out and the data were processed by K.M.A and M.M.S. The manuscript was written by T.C.L with supervision from P.A.R, except for §§2.3–2.6, which were written by K.M.A.

## Competing interests

We declare we have no competing interests.

## Funding

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

## Acknowledgements

The scans were performed at the University of New South Wales.

## Appendix A. Methods used to calculate responses

The model outlined in §2.1 is expressed in terms of a set of nonlinear partial differential equations with three spatial dimensions plus time. The main nonlinear term is due to the constitutive relation between pressure and blood mass density, which gives the term in equation (2.5) that depends on the blood mass density to the power of *β*, instead of just the blood mass density. This allows for nonlinear wave effects, and means that solution methods which depend on the equations being entirely linear, such as using Fourier transforms, as has been done previously by Aquino *et al.* [3], will not work. The three spatial dimensions present also mean that if the spatial resolution is increased, the computational costs for calculating the responses will increase with the cube of the number of points in any particular dimension. This means that the method for solving the equations should not require the spatial resolution to be too fine, or else it will take a prohibitively long time to calculate the response to a given drive.

The original way that model calculations were performed by Aquino *et al.* [3] was by linearizing the model, and then using Fourier transforms to turn the set of differential equations into a series of transfer functions in Fourier space that, combined with a standard fast Fourier transform algorithm, allowed the equations to be solved. The calculations performed this way predicted that there would be significant spreading of the BOLD signal for a localized neural drive, due to wave propagation, and this was confirmed in experiments [6].

In the nonlinear case, linear Fourier transform approaches are not appropriate. Instead, a finite differencing scheme is used for all the equations which relates the value of partial derivatives at a grid point to the values of quantities at that grid point and adjacent ones. For example, if the location in (*t*, *x*, *y*, *z*) space of a grid point is (*a*, *b*, *c*, *d*), then to calculate the second time derivative of , it uses
rmA1

where Δ*t* is the time step, and the *O* indicates that the error term is no larger than a constant multiplied by (Δ*t*)^{2}. Similarly, we can calculate the second spatial derivative in the *x* direction of using
rmA2

where Δ*x* is the grid spacing in the *x* dimension. The Courant–Friedrichs–Lewy (CFL) condition,
rmA3involves the speed *ν _{β}* at which waves propagate through the grid to the spatial (Δ

*x*, Δ

*y*, Δ

*z*) and temporal spacing (Δ

*t*) of those grid points. As long as this condition is satisfied, the equations can be solved by applying a nonlinear two-step Adams–Bashforth method [32], which has an error of

*O*(

*h*

^{2}) at each step, where

*h*is the largest of Δ

*t*, Δ

*x*, Δ

*y*and Δ

*z*. Continuing this for

*n*steps, where

*n*is proportional to 1/(Δ

*t*), will give an overall error in the method of

*O*(

*h*).

In our work, the initial conditions are set so that the blood mass density profile in the direction perpendicular to the cortical sheet matches the profile found by Lauwers *et al.* [14], and with the pressure and dHb content of the blood matching their equilibrium values, which gives an initial BOLD response of 0. The boundary conditions for the differential equations are open in the direction perpendicular to the sheet (because the blood mass density is able to vary at any point in the depth profile), but are closed on the other boundaries. However, the other boundaries are sufficiently far away from the effects the stimuli induces that any disturbance decreases to a negligible level by the time it reaches the edge.

Although the global error could be reduced by implementing a more complicated solution algorithm (using Runge–Kutta to increase the order of the error term), the error is already very small compared with the calculated value of the response (less than 1% of the size of the response) and so doing this would increase the accuracy only very slightly. In addition, there are much larger sources of error in the model for such small length and time scales, mainly from the smoothing mean field approximations to the actual physical system, so the increase in accuracy would not be useful anyway.

Using an explicit solving method (Adams–Bashforth) as opposed to an implicit method means that the stiffness of the differential equation does need to be taken into consideration. However, as the solutions were always physically reasonable and agreed with the results calculated using the linearized form of the equations, we find that divergence of the solution is not a problem for the time scales being used here.

In particular, the model equations are solved in the following order, in each case by replacing derivatives with the values at adjacent grid points as in equation (A 1) and rearranging to solve for the values that correspond to future times (in equation (A 1), this corresponds to the term):

(i) Equation (2.1) is solved to determine

*F*(the cerebral blood flow) from the neural drive*ζ*for the next time step.(ii) Equation (2.5) is solved to determine (the cerebral blood volume) and hence

*P*(the pressure, from equation (2.2)) from*F*for the next time step, taking care to satisfy the CFL condition while solving.(iii) Equation (2.6) is solved to determine

*Q*(the level of dHb) from and*P*for the next time step.(iv) Equation (2.7) is solved to determine

*Y*(the BOLD response) from and*Q*for the next time step.

In previous papers [6], the *z*-direction was averaged over, which explicitly assumed that the dynamics in that direction were separable from the rest. Although this is often a good approximation, it cannot always hold precisely, due to haemodynamic waves requiring time to propagate between the cortical layers. Using equation (A 1) to calculate the dynamics in the *z*-direction avoids this issue, so this is another benefit of implementing the new algorithm.

- Received July 25, 2016.
- Accepted November 17, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.