A modified damped Richardson–Lucy algorithm to reduce isotropic background effects in spherical deconvolution

Abstract

Spherical deconvolution methods have been applied to diffusion MRI to improve diffusion tensor tractography results in brain regions with multiple fibre crossing. Recent developments, such as the introduction of non-negative constraints on the solution, allow a more accurate estimation of fibre orientations by reducing instability effects due to noise robustness. Standard convolution methods do not, however, adequately model the effects of partial volume from isotropic tissue, such as gray matter, or cerebrospinal fluid, which may degrade spherical deconvolution results. Here we use a newly developed spherical deconvolution algorithm based on an adaptive regularization (damped version of the Richardson–Lucy algorithm) to reduce isotropic partial volume effects. Results from both simulated and in vivo datasets show that, compared to a standard non-negative constrained algorithm, the damped Richardson–Lucy algorithm reduces spurious fibre orientations and preserves angular resolution of the main fibre orientations. These findings suggest that, in some brain regions, non-negative constraints alone may not be sufficient to reduce spurious fibre orientations. Considering both the speed of processing and the scan time required, this new method has the potential for better characterizing white matter anatomy and the integrity of pathological tissue.

Introduction

In recent years diffusion MRI techniques have become an important research tool at the forefront of clinical neurosciences (Jones, 2008, Le Bihan et al., 2001). Diffusion tensor imaging (DTI) (Basser et al., 1994) and related tractography techniques (Basser et al., 2000, Behrens et al., 2003, Conturo et al., 1999, Mori et al., 1999, Parker et al., 2003) are the best known of these techniques. With these methods it is possible to study, in vivo, the microstructural organization of brain tissue and to extract several quantitative indices of anatomical integrity (Basser and Pierpaoli, 1996, Le Bihan, 2003). These models are not without limitations however. In particular the diffusion tensor model is unable to resolve multiple white matter orientations within voxels (Alexander, 2006, von dem Hagen and Henkelman, 2002). Alternative methods, such as diffusion spectrum imaging (DSI) (Wedeen et al., 2005), Q-ball imaging (Tuch, 2004), multiple compartment modelling (Behrens et al., 2007, Behrens et al., 2003, Jbabdi et al., 2007), Bayesian random effect modelling (King et al., 2009), Hybrid (Wu and Alexander, 2007), persistent angular structure MRI (PASMRI) (Jansons and Alexander, 2003) and spherical deconvolution (SD) (Anderson, 2005, Tournier et al., 2004) have been developed to address this problem.

Among these, the SD approach has the advantage of relatively short acquisition times, which are close to standard DTI clinical protocols, reduced computational times compared to some of the other methods, and the ability to resolve crossing fibres with a good angular resolution. SD is based on the assumption that the acquired diffusion signals from a single voxel can be modelled as a spherical convolution between the fibre orientation distribution (FOD) and the fibre response function that describes the common signal profile from the white matter (WM) bundles contained in the voxel (Tournier et al., 2004). Hence, the FOD could, in theory, provide a description that is closer to the physical fibre orientations than other methods that rely on the estimation of the water molecular displacement profile or the spin propagator (Alexander, 2005a, Tournier et al., 2007).

With the spherical deconvolution approach, highly ill-conditioned solutions are likely to occur due to the ill-posed nature of the deconvolution problem however. In the absence of regularization of the solution, even small changes in the acquired signal (e.g., MR noise) can lead to non-physical results (Bertero and Boccacci, 1998, Tournier et al., 2004). The introduction of non-negative constraints to the SD approach can improve the FOD estimation (Alexander, 2005a, Dell'Acqua et al., 2007, Jian and Vemuri, 2007, Tournier et al., 2007). Despite these recent developments, other sources of instability can still significantly affect the results in SD. The signal contribution from isotropic tissue, for example, is usually not included in spherical deconvolution models. Some authors have included an isotropic compartment in their signal model but these methods either require multiple b-value acquisitions (Jespersen et al., 2007), or apply stronger constraints on the solution, such as a finite number of fibre orientations (Behrens et al., 2007, Behrens et al., 2003) or known diffusivities (Kaden et al., 2008). No studies have, to date, investigated the direct effect of isotropic contributions on the FOD estimation.

To exclude voxels with isotropic partial volume (e.g., from grey matter (GM) or cerebrospinal fluid (CSF)), a threshold on fractional anisotropy maps (Pierpaoli and Basser, 1996), or similar maps, can been used to reject voxels with low anisotropy values. This approach is not a viable solution in those regions with low anisotropy due to highly complicated fibre orientations. Alternatively, Sakaie and Lowe have recently proposed the use of an optimised level of regularization of the estimated FOD at each voxel (Sakaie and Lowe, 2007).

In this study we combine a previously published Richardson–Lucy (RL) spherical deconvolution algorithm (Dell'Acqua et al., 2007) with an adaptive regularization technique to reduce spurious and non-physical fibre orientations in regions affected by partial volume whilst, at the same time, preserving angular resolution of the main fibre orientations. This new iterative algorithm can be described as a damped version of the RL–SD algorithm (dRL) as it uses the absolute dynamic range of the recovered fibre orientation to modulate the regularization of the solution. The results of simulations and in vivo data analysis are compared to assess the performance of the damped and the standard non-negative constrained RL algorithm using different levels of isotropic partial volume. The potential of the proposed method for application to patient studies is also assessed using scanning parameters that are usually applied in conventional DTI analysis (Jones et al., 2002).