Estimation of Laminar BOLD Activation Profiles using Deconvolution with a Physiological Point Spread Function

Background The specificity of gradient echo (GE)-BOLD laminar fMRI activation profiles is degraded by intracortical veins that drain blood from lower to upper cortical layers, propagating activation signal in the same direction. This work describes an approach to obtain layer specific profiles by deconvolving the measured profiles with a physiological Point Spread Function (PSF). New Method It is shown that the PSF can be characterised by a TE-dependent peak to tail (p2t) value that is independent of cortical depth and can be estimated by simulation. An experimental estimation of individual p2t values and the sensitivity of the deconvolved profiles to variations in p2t is obtained using laminar data measured with a multi-echo 3D-FLASH sequence. These profiles are echo time dependent, but the underlying neuronal response is the same, allowing a data-based estimation of the PSF. Results The deconvolved profiles are highly similar to the gold-standard obtained from extremely high resolution 3D-EPI data, for a range of p2t values of 5-9, which covers both the empirically determined value (7.1) and the value obtained by simulation (6.3). Comparison with Existing Method(s) Corrected profiles show a flatter shape across the cortex and a high level of similarity with the gold-standard, defined as a subset of profiles that are unaffected by intracortical veins. Conclusions We conclude that deconvolution is a robust approach for removing the effect of signal propagation through intracortical veins. This makes it possible to obtain profiles with high laminar specificity while benefitting from the higher sensitivity and efficiency of GE-BOLD sequences.


INTRODUCTION 42
Within the iso-cortex the type and density of brain cells varies following a laminar pattern. Following the 43 canonical model (Felleman and Van Essen, 1991) there are distinctive patterns of laminar connectivity, 44 corresponding to feedforward, feedback, and lateral interactions. Therefore, detecting neuronal activity at 45 the laminar level is of interest to understand the hierarchical relationship between different brain regions 46 and better understand brain function.

47
Although some initial attempts have been made using magnetoencephalography (Bonaiuto et al., 2018; 48 Troebinger et al., 2014), functional MRI (fMRI) is to date the most popular non-invasive in-vivo technique 49 to measure brain function at sub-millimetre resolution and obtain cortical depth-resolved activation signals

51
The most widely used method in fMRI measures the Gradient Echo (GE) based blood oxygenation level 52 dependent (BOLD) signal. The advantage of this method is that it has a higher sensitivity to activation and 53 a higher efficiency compared to the other MR methods available for fMRI. The disadvantage is its reduced 54 spatial specificity due to the contribution from large veins (Boxerman et al., 1995). In layer-specific fMRI the 55 spatial specificity of the measured BOLD profile is degraded by the contribution from emerging or 56 intracortical veins. These are the veins that are oriented perpendicular to the pial surface and drain blood 57 uni-directionally from lower to upper layers. Therefore, part of the activation signal in a lower layer will 58 propagate or leak into upper layers downstream. This phenomenon will be referred to as the "inter-laminar 59 leakage-problem" in this manuscript.

60
There are a number of MR methods that suffer less from large vein contamination. These were originally 61 developed for use in standard-resolution fMRI, but have been applied in laminar fMRI studies too. One of 62 these methods consists of combining the standard GE-BOLD method with a differential paradigm (Cheng  . The signal resulting from the subtraction of the condition of interest and the control condition is change in vascular blood volume. As the largest portion of blood volume change is located in the arterioles 77 and capillary bed, (Krieger et al., 2012) this technique is not weighted towards larger post-capillary vessels.

78
The alternative methods to GE-BOLD mentioned above present different specific technical difficulties of 79 their own (for example: sensitivity to motion, risk of exceeding energy deposition limits at high static 80 magnetic field strengths, inefficient tagging or inversion of the blood signal). Nonetheless, they all have in 81 common a lower sensitivity to activation and temporal resolution compared to the standard GE-BOLD based 82 fMRI approach.

83
This manuscript examines a method to deal with the inter-laminar leakage problem in the steady state. It

92
do not address the same problem. The deconvolution approach aims at removing the effect of draining 93 veins from the functional profiles whereas spatial GLM is a purely geometric approach that aims to reduce 94 partial volume effects, and is not necessarily related to functional acquisitions.

95
The deconvolution approach has previously been proposed  and 96 subsequently adopted by other users (Marquardt et al., 2018). It is attractive because although 97 deconvolution can be an unstable process, here it is equivalent to a subtraction of responses in lower layers 98 from the responses in higher layers, which is mathematically more robust. By use of a more sophisticated 99 modelling approach it can also be extended to the dynamic situation (Havlicek and Uludag, 2020). In the 100 current contribution the validity of the deconvolution approach is assessed on previously published profiles 101 acquired at 7T. First the approach is validated on the profiles obtained using an ultrahigh resolution 3D-EPI 102 acquisition approach (Fracasso et al., 2018). The profiles obtained in this study were grouped into strong, 103 middle and weak linear trends, whereby the latter were assumed to have a very small contribution from 104 intracortical veins. These latter profiles are considered as the gold-standard of leakage-free activation 105 profiles in this manuscript. Successful deconvolution of the full profile should hence yield a profile similar to 106 those with little contribution from intracortical veins. We then deconvolve the profiles at different TEs 107 obtained using a multi-echo 3D-FLASH sequence (Koopmans et al., 2011). In this dataset, BOLD profiles 108 are echo dependent, mostly due to the strong TE-dependency of the intravascular venous contribution, but 109 the underlying neuronal state is not. This feature is exploited to perform a data-based estimation of the PSF.

110
Lastly, the sensitivity of the deconvolution results to the choice of parameters used for the PSF is studied.

Mathematical description of the inter-laminar leakage 114
In the steady state, the measured activation profile can be understood as the weighted sum of the laminar 115 BOLD responses convolved with the corresponding layer specific physiological point spread function (PSF).

116
This can be expressed as Y = X*β + ε, where X is a matrix of layer specific PSFs, β is the underlying laminar 117 BOLD response, and ε is the error term. Specifically, the elements of the equation Y = X*β + ε are: where yi is the BOLD response measured at bin i, peaki is the peak value of the PSF measured at bin 119 (BOLD response directly related to the local neuronal response), spreadi,k is the spread of the PSF 120 measured at bin that results from the activation upstream in bin . β is the underlying laminar BOLD  2.
The peak and tail values might vary between subjects due to the inter-subject variability of the BOLD 158 response. However, given that the vascular architecture is similar between individuals, their ratio (the peak 159 to tail ratio, p2t) can be expected to show less variance. There are two ways of determining the elements 160 of matrix X in Equation 2. One of the options, which we will refer to as normalised, is to set peak = y1 and 161 tail = y1/p2t, in which case the underlying profile of the BOLD response, β, will be relative to the response 162 in the lowest bin. The other option is to write the matrix as all ones on the diagonal and 1/p2t in the lower 163 diagonal. In this second approach, the laminar BOLD profiles will not be normalised and differences in the 164 average magnitude of the profiles between subjects would reflect differences in the average BOLD 165 response between them, and variations in the profiles would represent both the underlying neuronal activity 166 and differences in the underlying vascular density.

167
The consequences for the laminar profiles of: using an exact PSF for each layer; normalising this exact 168 PSF to the response at the lowest bin; or replacing the normalised, layer-specific PSFs with a single 169 normalised p2t value that is constant across the cortex, are described in the Supporting Information (part

186
The p2t values presented in the previous paragraph assume that the cortex has been divided into 10 bins.

187
If a different number of bins is used, then the p2t has to be adjusted according to the following formula: where n is the number of bins across the cortex in the acquired data and nmodel is the number of bins used

197
PSF in matrix X will vary with TE (as described below) but the pattern followed by the underlying laminar 198 response, β, will be TE-independent. Therefore, when activation profiles at different echo times are 199 combined, the system will not be underdetermined and both the TE-dependent peak and tail magnitudes

200
(and hence the p2t) and the TE-independent lamina BOLD responses can be estimated. The multi-echo 201 3D-FLASH sequence used had 10 echo times ranging from TE= 5.7 ms to TE= 56.1 ms (see next 202 subsection for more details on these data). In order to better compare profiles acquired at different TEs it 203 was found necessary to also model the signal contribution from the pial leading to a slight extension of 204 Equation 2. This is particularly important for short TEs where the intravascular contribution at the pial 205 surface can be large. Two parameters, the pial response, , and the partial volume factor, , are 206 introduced to account for partial volume effects with the pial surface signal as follows: Equation 4 where is the fraction of pial partial volume and can spread up to two bins (equivalent to one voxel) below

228
 Peak values: randomly picked for each TE from a uniform distribution between 0 and 10.

229
 Tail values: the peak value divided by a plausible peak to tail ratio (p2t). This peak to tail ratio was 230 randomly chosen from a uniform distribution between 4 and 8. This range was chosen based on 231 the results presented in Figure 5, and because the p2t estimated in the model was 6.3.

232
 Pial values: randomly picked for each TE from a uniform distribution between 0 and 10.

233
As no constraints on the profiles of the βs were set, very large and small values of p2t could also represent 234 numerical solutions for the system, although they are not physiologically plausible. That is why a loose used for the deconvolution are given (see Table 1 for a summary).

250
The ultrahigh-resolution dataset was acquired at 7T using the 3D-EPI sequence and a dedicated coil array   If the similarity metric is 1, the two profiles are identical or scaled versions of one another. If 0, they are 288 orthogonal to each other.

289
In the present article we consider the weak linear trend profiles to represent a gold standard for the

301
Deconvolution with the p2t on GE-BOLD activation profiles

302
The deconvolution results for the ultrahigh resolution dataset are shown in Figure 2. Figure 2a shows the 303 weighted average cortical profile for each subject (one acquired twice), obtained using Equation 6, the 304 average across all subjects is shown in black. Figure 2b shows the corresponding deconvolved profiles.

305
The ultrahigh resolution dataset was sampled with 16 bins across the cortex, therefore, before performing

336
In the multi-echo dataset, measured activation profiles and p2t values are expected to be TE-dependent, 337 but the underlying neuronal activity is not. This feature is exploited to obtain a data-based estimation of the 338 p2t by combining the data for all echo times and solving the system in Equation 5 following a minimisation 339 approach. Figure 7 shows the results of the p2t estimation for TE=27.6 ms, close to the T2* of grey matter 340 at 7 T. Most of the p2t values obtained fall between 5 and 9. The results for Subject 1 represent an outlier 341 that can be explained by the unusual shape of the profiles for this subject at longer TEs. When these longer 342 echo times, with a poorer SNR, were not considered, then Subject 1 was no longer an outlier and the 343 distribution of the estimated p2ts was concentrated more strongly in the range 5-9 (SI: Figure 5).

363
The deconvolved profiles in the multi-echo dataset tend to show a slight increase just above the middle of 364 the cortex, where microvascular density shows its maximum (Duvernoy et al., 1981)    will be accompanied by differences in the intracortical veins that drain them. The form of the PSF is not 379 expected to vary and the fact that the parameter that describes the PSF is used as a ratio in this approach 380 will minimise differences related to differences in the vascular density between regions.

381
The results in this manuscript show that a peak to tail ratio ~6.3 for TE= 27.6 ms at 7T for 10 sampling 382 points across the primary visual cortex adequately describes the biophysical properties of the leakage effect.

383
If the same approach is applied to other regions then the p2t value will need to be adjusted to accommodate

394
General requirements to apply the deconvolution approach 395 Most importantly good segmentation and co-registration are paramount to applying the approach. If the 396 cortical sampling deviates too much from the model due to segmentation problems, the approach will fail.

397
Once the cortical boundaries have been defined it is also important that the bins obtained are the same 398 size. The model used to estimate the p2t,  implicitly relies on the equivolume 399 sampling approach (Waehnert et al., 2014). Hence, this method should be favoured over the equidistant 400 approach, although the differences between the two have been found to be negligible for voxels with a side > 401 0.6 mm in the direction perpendicular to the cortical surface (Kemper et al., 2017).

402
As previously mentioned Equation 2 assumes that the response of the cortical vasculature to the stimulus 403 is in a steady state and does not apply to transient cases, this will be the case for stimuli longer than 1s.

411
The multi-echo dataset was acquired using the 3D-FLASH sequence. Therefore, the BW per pixel in the 412 PE direction is infinite whereas the BW/pixel in the RO direction is 240Hz/pixel. These values are above 413 the line broadening induced by the acquisition train and hence the results presented here will not be affected.

414
The high-resolution dataset was acquired using 3D-EPI and the BW/pixel in the phase encoding directions 415 was 21Hz/pixel (Fracasso et al., 2018). This is smaller than the line broadening expected for venous blood.

416
However, the TE= 27ms was matched to T2*GM. At this echo time, the intravascular contribution to the 417 BOLD signal is more than 10 times smaller than the extravascular contribution (Cheng et al., 2015). Hence,

418
it can be assumed that the smearing of the profiles due to Lorentzian line-broadening is negligible.

419
Regardless of whether deconvolution is used or not, studies focusing on laminar activation or that require 420 high spatial accuracy should control that the Lorentzian broadening is confined to a voxel when setting the 421 acquisition parameters.

424
The deconvolution approach has been applied to GE-BOLD profiles obtained at 7T acquired using similar