Patent Application
20230160987
This invention relates to magnetic resonance imaging, and more particularly to quantitatively mapping material intrinsic physical properties using signals collected in magnetic resonance imaging.
Quantitative susceptibility mapping (QSM) in magnetic resonance imaging (MRI) has received increasing clinical and scientific interest. QSM has shown promise in characterizing and quantifying chemical compositions such as iron, calcium, and contrast agents including gadolinium and superparamagnetic iron oxide (SPIO) nanoparticles. Tissue composition of these compounds may be altered in various neurological disorders, such as Parkinson's disease, Alzheimer's disease, stroke, multiple sclerosis, hemochromatosis, and tumor as well as other disease throughout the body. QSM is able to unravel novel information related to the magnetic susceptibility (or simply referred as susceptibility), a physical property of underlying tissue. Due to the prevalence of iron and calcium in living organisms, their active involvement in vital cellular functions, their important roles in the musculoskeletal system, QSM is in general very useful to study the molecular biology of iron/calcium by tracking iron in the circulation system and the metabolic activities by using iron and calcium as surrogate marks. Therefore, accurate mapping of the magnetic susceptibility induced by iron, calcium and contrast agents will provide tremendous benefit to clinical researchers to probe the structure and functions of human body, and to clinicians to better diagnose various diseases and provide relevant treatment.
Implementations of systems and methods for collecting and processing MRI signals of a subject, and reconstructing maps of physical properties intrinsic to the subject (e.g., magnetic susceptibility) are described below. In some implementations, MR signal data corresponding to a subject can be transformed into susceptibility-based images that quantitatively depict the structure and/or composition and/or function of the subject. Using this quantitative susceptibility map, one or more susceptibility-based images of the subject can be generated and displayed to a user. The user can then use these images for diagnostic or experimental purposes, for example to investigate the structure and/or composition and/or function of the subject, and/or to diagnose and/or treat various conditions or diseases based, at least in part, on the images. As one or more of the described implementations can result in susceptibility-based images having higher quality and/or accuracy compared to other susceptibility mapping techniques, at least some of these implementations can be used to improve a user's understanding of a subject's structure and/or composition and/or function, and can be used to improve the accuracy of any resulting medical diagnoses or experimental analyses.
Some of the implementations described below can be used to evaluate the magnetic dipole inversion while allowing incorporation of preconditioning, anisotropic weighting, primal dual solver, chemical composition based fat compensation, and enforcement of regional known near-constant susceptibility is disclosed to improve susceptibility map quality including suppression of artificial signal variation in homogenous regions and to provide absolute susceptibility values, which are pressing concerns to ensure quantification accuracy, experimental repeatability and reproducibility.
In general, one aspect of the invention disclosed here enables an absolute quantification of tissue magnetic susceptibility by referencing to a region of tissue with known homogenous susceptibility, such as zero susceptibility for cerebrospinal fluid (CSF) in the brain and fully oxygenated aorta in the body. This aspect of invention solves the practical problem of zero-reference in longitudinal and cross-center studies and in studies requiring absolute quantifications. The invention disclosed here also enables reduction of shadow artifacts caused by field data that are not compatible with the dipole field model, in addition to reducing streaking artifacts caused by granular noise in data.
In general, another aspect of the invention disclosed here enables accurate tissue structure matching using an orientation dependent L1 norm and a primal-dual computation of its derivative in numerical searching for the best match. The disclosed preconditioned total field inversion method bypasses the background field removal step, the robustness of which may be affected by phase unwrapping, skull stripping, and boundary conditions.
In general, another aspect of the invention disclosed here enables accurate imaging methods to quantify liver iron by correcting for effects from fat, fibrosis and edema that confound traditional T2 and T2* methods, to quantify cerebral metabolic rate of oxygen consumption by incorporating cerebral blood flow information, to quantify bone mineralization and other calcifications in the body, and to quantify chamber blood oxygenation level and intramyocardial hemorrhage in the heart.
Implementations of these aspects may include one or more of the following features.
In some implementations, the cerebrospinal fluid in the ventricles in the brain is segmented based on low R2* values by applying R2* thresholding and connectivity.
In some implementations, the prior information regarding the magnetic susceptibility distribution is determined based on an L1 norm and structural information of the object estimated from acquired images of the object.
In some implementations, the object comprises a cortex, a putamen, a globus pallidus, a red nucleus, a substantia nigra, or a subthalamic nucleus of a brain, or a liver in an abdomen, and generating one or more images of the object comprises generating one or more images depicting the cortex, the putamen, the globus pallidus, the red nucleus, the substantia nigra, or the subthalamic nucleus of the brain, or the liver in the abdomen.
In some implementations, the object comprises at least one of a multiple sclerosis lesion, a cancerous lesion, a calcified lesion, or a hemorrhage, and generating one or more images of the object comprises generating one or more images depicting at least one of the multiple sclerosis lesion, the cancerous lesion, the calcified lesion, or the hemorrhage.
In some implementations, the operations further comprise of quantifying, based on the one or more images, a distribution of contrast agents introduced into the object in the course of a contrast enhanced magnetic resonance imaging procedure.
In some implementations, tissue chemical composition is fat and water, and the magnetic field, the water fraction and the fat fraction in a voxel are computed iteratively from the complex MRI data using an appropriate initialization.
The details of one or more embodiments are set forth in the accompanying drawings and the description below. Other features and advantages will be apparent from the description and drawings, and from the claims.
Like reference symbols in the various drawings indicate like elements.
Implementations of systems and methods for collecting and processing MRI signals of a subject, and reconstructing maps of physical properties intrinsic to the subject (e.g., magnetic susceptibility) are described below. Some of the disclosed implementations overcome limitations in quantitative susceptibility mapping (QSM), where the computation speed and accuracy and the Bayesian prior information affect QSM image quality and practical usage. In some implementations, to reconstruct maps of a subject's magnetic susceptibility, the iterative computation process of finding the susceptibility distribution that optimally fit available field data and tissue structure prior information is accelerated using preconditioning and primal dual optimization method and is computed more accurately using primal dual optimization method. In some implementations, the tissue structure prior information in a Bayesian method for estimating susceptibility is improved using anisotropic edge weighting and using additional information of the near-zero susceptibility of the cerebrospinal fluid (CSF) in the lateral ventricles that can be automatically segmented from the brain image. Implementations of these improvements have made QSM technology robust and accurate, extending QSM usage from the brain to including the skull, bone, heart and liver, and to quantification of deoxyheme associated with aerobic energy metabolism.
QSM is based on minimizing a cost function that fits the measured tissue magnetic field b(r) (in unit B0) to the dipole field
Approximating noise in b as Gaussian, the data fidelity term is
w(r)|d*χ(r)−b(r)|2dr≡∥d*χ−b∥w2 [1]
This integral form b=d*χ allows numerical estimation of the susceptibility solution, typically under a regularization associated with a prior information. As known from physics of Maxwell equation, the integral form may represent a solution to the governing partial differential equation (PDE). For mathematical analysis of susceptibility solution, the equivalent PDE that can be derived directly from Maxwell's equation with Lorentz sphere correction offers useful insights:
Eq.2 is a wave equation with respect to susceptibility χ (z-axis is time) and a Laplacian equation with respect to field b. In the continuous space without any error with b(r)=d*χ(r), B(k)=D(k)X(k), a precise solution of no streaking can be obtained: χ=((1−η∈)(B/D)). η∈=1 if D(k)<∈/2 otherwise 0. Any field data that cannot be fitted to the dipole field is referred as the incompatible field data v(r), which may come from various causes including noise, digitization error, and anisotropy source. This incompatible field causes artifacts according to the wave propagation solution:
is the wave propagator having large values at magic-angle cones above and below the incompatible source
The dipole incompatible part, including noise, discretization error, and anisotropy sources, generates artifacts defined by the wave propagator with the B0 field direction as the time axis. Granular noise causes streaking and contiguous error causes both streaking and shadow artifacts.
The invention described here makes use of the above mathematical fact to improve QSM speed and accuracy. Streaking artifacts are characterized by edges along the magic angle (54.7°) in k-space or the complimentary magic angle in image space, which are almost completely distinct from tissue edges. Accordingly, they can be minimized during numerical optimization by a regularization term based upon tissue structure information to identify solution of no streaking, which is the morphology enable dipole inversion (MEDI) approach. Tissue structure information is numerically evaluated using edge gradient and L1 norm. This calculation was previously performed by approximating the derivative of L1 norm as a continuous function and by ignoring the edge gradient direction. This approximation can be improved in terms of both accuracy and speed using the primal dual optimization method and using anisotropic weighting according to edge gradient direction.
Shadow artifacts have low spatial frequency contents, which can be minimized accordingly by incorporating cost function terms of low spatial frequencies. Additionally, penalty terms can be formulated using specific tissue structure and shadow structure information. One specific tissue structure is that CSF in the ventricles of the brain is almost pure water with very little cellular contents. Therefore, ventricular susceptibility map should be nearly uniform, and any deviation from uniformity should be regarded as shadow artifacts to be penalized during numerical optimization by incorporating a regularization term. One specific shadow structure information is that shadow generated by smoothly distributed error source over a large region, such as white matter anisotropic component, which is strong above and below the region, spreads cross the whole region, and consists of components of low spatial frequencies. Therefore, the shadow artifacts can be penalized during numerical optimization using a k2 weighting (Laplacian operation) as a preconditioner to avoid low spatial frequency in Krylov space, or using a regularization term imposing costs on low spatial frequencies with corresponding long spatial wavelengths such as object size and half object size and a cost function associated with an estimated white matter anisotropic contribution. In another implementation, the white matter anisotropy effect modeled according to a brain atlas as a dipole-incompatible part is tuned such that the CSF susceptibility in the segmented ventricle region is uniform, which at the same time minimizing shadow artifacts in the brain QSM.
The iteration for finding the optimal susceptibility distribution according to available field data and tissue prior information can be hugely accelerated using preconditioning. This is particularly beneficial for searching in a large range of possible susceptibility values. As a result, the improvement in QSM system and method described in this invention enables fast, accurate and robust QSM, which improves applications in the brain and enables applications in other organs outside the brain.
Robust QSM improves upon many applications, including brain QSM with reduced shadow and streaking artifacts and with uniform ventricle CSF as an automatic and consistent zero reference, whole brain QSM without skull stripping and bone mineralization where bone diamagnetic susceptibility being proportional to CT attenuation coefficients can be used for positron emission tomography attenuation correction, liver iron content mapping without R2/R2* confounding errors, mapping cerebral metabolic rate of oxygen consumption, mapping susceptibility sources in the heart including oxygenation level and intramyocardial hemorrhage, and mapping susceptibility sources in vessel wall plaque differentiating intraplaque hemorrhage from plaque calcification.
In this section, we describe a QSM improvement that uses preconditioning to speed up the convergence of QSM, allowing integration of background field removal and tissue field inversion.
Quantitative susceptibility mapping (QSM) extracts the spatial distribution of tissue magnetic susceptibility from the gradient echo signal phase. Current QSM methods consist of two steps: i) background field removal to determine the local field generated by tissue, and ii) inversion from the local field to the tissue susceptibility. This allows a fairly robust mapping of the central brain regions, particularly for iron deposition in the mid brain nuclei. However, several technical challenges remain.
A major challenge is imprecise separation of background and tissue fields caused by the assumptions made in background field removal methods. This is particularly problematic near the brain boundary where a large tissue-air susceptibility difference is present. To avoid the separate fitting of background and local field, Laplacian-based QSM methods have been proposed based on the partial differential formulation of the forward signal equation. The Laplacian operation implicitly eliminates the background field. However, the practical implementation of the Laplacian requires a trade-off between robustness to error amplification and the integrity of the visualized cortical brain tissue. The necessary erosion of the brain mask may prevent visualization of certain structures at the brain boundary, such as the superior sagittal sinus.
The presence of a large susceptibility dynamic range within the region of interest (ROI) is similarly a challenge in QSM, often leading to streaking artifacts. For example, the susceptibility difference between intracerebral hemorrhage (ICH) and surrounding tissue can exceed 1.6 ppm. Using a nonlinear QSM model of ICH signal can reduce these streaking artifacts, but does not eliminate them. Recent work has focused on this challenge by separating the fitting processes for sources of strong (such as ICH) and weak (normal brain) susceptibilities, hence preventing ICH-related artifact from permeating into the normal brain. However, these methods require carefully choosing the regularization parameter or threshold for detecting ICH to minimize artifacts in the subsequent inversion for weak susceptibilities.
In this work, we address both challenges using a generalized inverse problem that is suitable for large dynamic ranges in susceptibility. We propose a preconditioned QSM to perform total field inversion (TFI), strongly reducing the error propagation associated with imprecise background field removal, and suppressing streaking artifacts in intracerebral hemorrhage QSM.
For QSM, the total magnetic field is conventionally decomposed into two components: the local field ƒL defined as the magnetic field generated by the susceptibility χL inside a given region of interest M (local susceptibility), and the background field ƒB defined as the magnetic field generated by the susceptibility χB outside M (background susceptibility):
ƒ=ƒB+ƒL=d*(χB+χL) [4]
Here * is the convolution operator, and d is the field generated by a unit dipole with Lorentz sphere correction. The components ƒB and ƒL may be estimated separately: estimation of χB or ƒB is referred to as background field removal. A variety of background field removal methods have been proposed, such as high-pass filtering (HPF), projection onto dipole fields (PDF) or Laplacian based methods. χL is then obtained from the local field ƒ−ƒB, typically using regularization. This step is referred to as local field inversion (LFI). Errors in the background field propagate into the subsequent local field, ultimately leading to errors in the final susceptibility. Although recent work allows updating the background field during local field inversion, it requires a prior susceptibility atlas and co-registration to estimate the background susceptibility first.
Here, we propose to estimate χB and χL jointly using a total field inversion (TFI):
thus using the same formulation as the traditional QSM inversion problem. Here χ=χL+χB represents the total susceptibility in the entire image volume. The data weighting w can be derived from the magnitude images by calculating the error propagation from the MRI data into the total field. A weighted total variation term is used to suppress the streaking artifacts induced by the presence of zeroes in the dipole kernel d.
The use of iterative optimization algorithms, such as conjugate gradient (CG), in solving Eq. 5 can lead to slow convergence, as illustrated in a healthy brain in n and covariance matrix Γ∈
n×n, then a right-hand preconditioner P (see below) that approximates the covariance matrix Γ, or PHP≈Γ, will increase convergence speed. For TFI, a binary diagonal preconditioner P is constructed as follows:
where i denotes the voxel, and M the tissue ROI. It is designed such that the difference in the matrix PHP between a voxel outside and inside M is approximately equal to the difference in susceptibility between the local and the background region, which includes bone and air. We further modify this P to account for strong susceptibility within M (e.g. intracerebral hemorrhage) by thresholding the R2* map: MR2*:=R2*<R2th*, assuming that voxels with high R2* have strong susceptibility. Thus, P is defined as:
The preconditioned TFI problem then becomes:
The final susceptibility is computed as χ*=Py*. Eq. 8 is solved using a Gauss-Newton method as follows:
Where F is the Fourier transform,
is the k-space dipole kernel and G the finite difference operator. At the n-th Newton iteration, Φ is linearized at the current solution yn:
Using weak derivatives, the update dyn is found by solving:
which is solved using the conjugate gradient (CG) method. Note that we use MGGPyn≈√{square root over (|MGGPyn|2+∈)} where ∈=10−6 to avoid division by zero. The stopping criteria are: maximum number of iterations=100 or relative residual <0.01. The Gauss-Newton method is terminated when the relative update ∥dyn∥2/∥yn∥2 is smaller than 0.01.
A number of recently proposed methods similarly circumvent background field removal to estimate susceptibility directly from the total field. These methods are based on the partial differential formulation of the signal equation (Eq. 4):
Here, Δ denotes the Laplacian. Notice that the contribution from background sources χB disappears after applying the Laplacian. Two examples of these methods are:
where MG is a binary mask derived from the magnitude image. It can be efficiently optimized using a preconditioned conjugate gradient solver.
Laplacian-based methods enable QSM reconstruction within a single step. However, the practical implementation of Laplacian necessitates the erosion of the ROI M. The amount of erosion depends on the width of the kernel used for computing the Laplacian.
For developing and evaluating the proposed total field inversion (TFI) with preconditioning, we performed simulations, phantom imaging, in vivo imaging of the brain of healthy subjects using COSMOS as a reference and of patients suffering from hemorrhages. Finally, we explored whole brain QSM. The optimal PB in Eq. 7 was empirically determined to ensure that, for a given number of CG iterations, the Gauss-Newton solver arrived at a solution with minimal error with respect to a reference susceptibility. Then this PB value was used for other similar datasets. In this work, we used R2th*=30 s−1. The performance of the proposed TFI method in phantom and in vivo was compared with SSQSM and Differential QSM. In this work, the Spherical Mean Value (SMV) with variable kernel size was used to implement the Laplacian in SSQSM and Differential QSM, where the kernel diameter was varied from d to 3 voxels, and the ROI mask was eroded accordingly for these two methods. The kernel diameter d at the most central ROI was optimized as described below.
Accuracy of PDF+LFI and TFI. We constructed a numerical brain phantom of size 80×80×64. Brain tissue susceptibilities were set to 0 ppm, except for a single point susceptibility source of 0.1 ppm. Background susceptibilities were set to 0, and the total field was computed with the forward model ƒ=Md*χ. Both LFI and TFI without preconditioning (PB=1) were used, and the estimation for the point source χS was compared to the truth χST (0.1 ppm). For sequential inversion, both PDF and LBV were used for background field removal and MEDI for local field inversion. This process was repeated by moving the point source across the entire ROI, which generated an error map showing the spatial variation of the estimation error |χS−χST|/χST for each method. The regularization parameter λ was set to 10−3 for PDF+LFI, LBV+LFI and TFI.
Effect of the preconditioner P in TFI. A numerical brain phantom of size 256×256×98 was constructed with known susceptibilities simulating different brain tissues: white matter (WM) −0.046 ppm, gray matter (GM) 0.053 ppm, globus pallidus (GP) 0.19 ppm, and cerebrospinal fluid (CSF) 0. Background susceptibilities were set to 9 ppm to simulate air outside the ROI. The total field was computed from the true susceptibility map χT (
To examine the accuracy of our TFI method, a phantom with 5 gadolinium solution filled balloons were embedded in agarose (1%), with expected magnetic susceptibilities χref of 0.049, 0.098, 0.197, 0.394 and 0.788 ppm, respectively. The measured susceptibilities were referenced to the agarose background. The phantom was scanned at 3T (GE, Waukesha, Wis.) using a multi-echo gradient echo sequence with monopolar Cartesian readout (anterior/posterior) without flow compensation. Imaging parameters were: FA=15, FOV=13 cm, TE1=3.6 ms, TR=71 ms, #TE=8, ΔTE=3.5 ms, acquisition matrix=130×130×116, voxel size=1×1×1 mm3, BW=±31.25 kHz and a total scan time of 13 min.
Next, PDF+LFI, TFI, SSQSM and Differential QSM were performed on the same phantom data. Nonlinear field map estimation followed by graph-cut-based phase unwrapping was used to generate the total field. We optimized PB in Eq. 7, by repeating the TFI methods for a range of PB values and selecting the value that produced the smallest error between the mean estimated susceptibility χ for each balloon and its corresponding reference susceptibilities χref: Errphantom=Σi=15 ((χi−χrefi)/χrefi) at the end of the 300th CG iteration. The regularization parameter λ was determined by minimizing Errphantom for PDF+LFI, which was then also used for TFI. For SSQSM and Differential QSM, λ and the maximum kernel diameter d for the variable radius SMV were jointly determined by minimizing Errphantom
The brains of 5 healthy subjects were imaged using multi-echo GRE at 3T (GE, Waukesha, Wis.) with monopolar readout without compensation. All studies in this work were conducted with the approval of our institutional review board. Imaging parameters were: FA=15, FOV=24 cm, TE1=3.5 ms, TR=40 ms, #TE=6, ΔTE=3.6 ms, acquisition matrix=240×240×46, voxel size=1×1×3 mm3, BW=±62.5 kHz and a total scan time of 9 min. For each subject, 3 orientations were imaged, and the reference brain QSM χLref was reconstructed using COSMOS. For PDF+LFI, TFI, SSQSM and Differential QSM, only one orientation was used. The ROI was obtained automatically using BET. Data from one subject was used to optimize PB by minimizing the RMSE Errinvivo=∥χL−χLref∥2/√{square root over (n)} between the estimated susceptibility χL and the reference χLref at the end of the 300th CG iteration (n is the number of voxels inside the ROI. The obtained weight PB was then used for the remaining 4 subjects. The CSF was chosen as the reference for the in vivo susceptibility value in this work. Susceptibilities were measured within ROIs of the GP, putamen (PT), thalamus (TH) and caudate nucleus (CN) for quantitative analysis. The relative difference (χGP−χGPCOSMOS)/χGPCOSMOS was calculated for each subject, with χGP the mean susceptibility measurement of GP for each subject. The regularization parameter λ was chosen by minimizing Errinvivo for PDF+LFI, while TFI used the same λ. The maximum kernel diameter d and the regularization parameter λ were jointly determined by minimizing Errinvivo for SSQSM and Differential QSM.
We acquired human data in 18 patients, each with intracerebral hemorrhage (ICH), brain tumors that including cancerous lesions, and calcified lesions, using multi-echo GRE at 3T (GE, Waukesha, Wis.) with monopolar readout without flow compensation. Imaging parameters were: FA=15, FOV=24 cm, TE1=5 ms, TR=45 ms, #TE=8, ΔTE=5 ms, acquisition matrix=256×256×28˜52, voxel size=1×1×2.8 mm3, BW=±31.25 kHz and parallel imaging factor R=2. Total scan time was proportional to the number of slices (about 10 slices/min). For each case, we applied PDF+LFI, TFI, SSQSM and Differential QSM to estimate the QSM of brain with ICH. For TFI, Eq. 4 was used and PB was determined from the previous in vivo COSMOS experiment. The R2* threshold was R2th*=30 s−1, which, in preliminary studies, was empirically determined to effectively distinguish the hemorrhage site from surrounding brain tissue. The regularization parameter λ was chosen using L-curve analysis for PDF+LFI, SSQSM and Differential QSM, while TFI used the same λ as in PDF+LFI. The same kernel diameters d as in the COSMOS experiment were used here for SSQSM and Differential QSM, respectively. To quantify the shadow artifact around ICH, the mean susceptibility within a 5 mm wide layer surrounding each ICH was computed for PDF+LFI and TFI. The reduction in standard deviation (SD) in non-ICH regions was also computed as a measure for ICH-related artifact reduction. It was defined as R=(SD(χPDF|Mnon-ICH)−SD(χTFI|Mnon-ICH))/SD(χPDF|Mnon-ICH), where χ|Mnon-ICH denotes the susceptibilities in the non-ICH region.
Since the proposed method no longer separates background field removal and local field inversion, it is possible to generate a susceptibility map for the entire head by including tissues outside the brain (e.g. scalp, muscles and oral soft tissues) into the ROI M in Eq. 8. We acquired one healthy human brain data at 3T (GE, Waukesha, Wis.) using a multi-echo gradient echo sequence with monopolar 3D radial readout for large spatial coverage. Flow compensation was off. Imaging parameters were: FA=15, FOV=26 cm, TE1=1 ms, TR=34 ms, #TE=9, ΔTE=3 ms, acquisition matrix=256×256×256, voxel size=1×1×1 mm3, BW=±62.5 kHz, number of radial projections=30000, scan time of 15 min, and reconstruction using regridding. TFI was used with PB determined from the previous in vivo COSMOS experiment. The ROI M in Eq. 8 was determined by thresholding the magnitude image I:M:=I>0.1 max(I). The total field was estimated using a graph cut based concurrent phase unwrapping and chemical shift determination. As phase wrapping is an integer operation in unit of 2π and graph cut technique is best for handling such discrete integration operation. Water fat separation problem can be estimated in two steps: 1) generate an initial estimate of phase unwrapping, susceptibility field and chemical shift by simplifying a voxel composition as either water or fat, 2) input the initial estimates into a full water-fat signal model to fine-tune susceptibility and water-fat fraction through numerical optimization. As water-fat is not a convex problem, the first step based on graph cut provides a robust initialization for water-fat numerical optimization procedure to converge to a realistic solution.
Differential QSM was also applied using this ROI. For comparison, QSM of only the brain was also reconstructed with TFI using the mask obtained from BET. The regularization parameter λ was chosen using L-curve analysis for brain-only TFI, Differential QSM and whole head TFI.
All computations in this work were performed in MATLAB on a desktop computer with a 6-core CPU (Intel Core i7) and 32 GB memory.
Accuracy of PDF+LFI and TFI. In
Effect of preconditioner P in TFI.
The optimized regularization parameter λ was 2×10−2 for both PDF+LFI and TFI, 5×10−1 for SSQSM and 5×10−3 for Differential QSM. The optimized kernel diameter d was 13 for SSQSM and 3 for Differential QSM. The weight PB=10 was found to be optimal here (
In Vivo Imaging: Brain with Hemorrhage
All hemorrhage patient brain images were reconstructed using PB=30 for TFI. The average reconstruction time was 328 seconds for PDF+LFI, 325 seconds for TFI, 50 seconds for SSQSM and 344 seconds for Differential QSM. The regularization parameter λ was 10−3 for both PDF+LFI and TFI, 2×10−3 for SSQSM and 5×10−4 for Differential QSM. Maximum kernel diameter d was 17 for SSQSM and 3 for Differential QSM.
The whole head QSM was compared with the magnitude/phase image and brain QSM (
Our data demonstrate that total field inversion (TFI) for QSM eliminates the need for separate background field removal and local field inversion (LFI), and that preconditioning can accelerate TFI convergence. Compared to the traditional PDF+LFI approach, TFI provides more robust QSM in regions near high susceptibility regions such as those containing air or hemosiderin present in hemorrhages, within similar reconstruction time. It is also demonstrated that TFI is able to generate the whole head QSM without the need for brain extraction.
The sequential background field removal and LFI process the same data using the same Maxwell equations but require an assumption or regularization) to differentiate the background field and the local field. For the background field removal exemplified here using PDF, it is assumed that the Hilbert space L spanned by all possible local unit dipole fields ƒdL is orthogonal to the Hilbert space B spanned by all possible background unit dipole fields ƒdB. However, this orthogonality assumption is not valid when the local unit dipole is close to the ROI boundary, and will cause error in PDF. This error or similar error introduced by any other background field removal method propagates into the subsequent LFI and produces an inaccurate local susceptibility estimation. However, there is no need to break the problem of fitting MRI data with the Maxwell equations into two separate background field removal and LFI problems. Our proposed TFI eliminates this separation and avoids the associated error propagation. Improvement of TFI over PDF+LFI is shown in Simulation (
Preconditioning is necessary to achieve practical performance with TFI.
Since our preconditioner is designed for a large dynamic range of susceptibilities, it straightforwardly applies to intracerebral hemorrhage (ICH). As shown in
Laplacian-based methods, SSQSM and Differential QSM use the harmonic property of the background field ƒB (ΔƒB=0 within the ROI) to eliminate both phase unwrapping and background field removal, enabling local susceptibility estimation in a single step. SSQSM further speeds up by omitting the SNR weighting and using L2 norm for regularization. However, Laplacian-based methods suffer from brain erosion: The Laplacian is implemented using the finite difference operator or the spherical kernel operator, both requiring the ROI mask to be eroded. As a consequence, erosion of superior sagittal sinus can be seen in SSQSM and Differential QSM (
The capability of quantitatively mapping both diamagnetic and paramagnetic materials as demonstrated here for the whole head QSM can be applied to atherosclerotic imaging to resolve and define calcification (diamagnetic with negative susceptibility) and intraplaque hemorrhage (paramagnetic with positive susceptibility) in plaques found in carotid arteries, coronary arteries, aorta and other vessel walls. The multi-echo gradient echo 3D sequence can be modified to include EEG gating to reduce effects from pulsatile changes over the cardiac cycle, and the flow compensation can be included to minimize flow generated phase. The QSM of plaque calcification and intraplaque hemorrhage components would solve the problem of confounding calcification and hemorrhage, both of which appear hypointense in magnitude images, allowing measurement of intraplaque hemorrhage that is a major determinant of plaque vulnerability.
The whole head QSM provides skull bone calcification spatial distribution quantification, which has been shown to be linearly proportional to the attenuation coefficients measured on CT. This skull bone and other intracranial calcifications measured on QSM can be used to correct for ultrasound attenuation and scattering in magnetic resonance guided focused ultrasound (MRgFUS), and to correct for gamma attenuation MR/PET that simultaneously forms MRI and positron emission tomography (PET) images. Both MRgFUS and MR/PET requires correction based on quantitative bone spatial distribution, which are not available from standard MRI and may require an additional CT imaging session with costs and radiation. Bone defined on QSM can eliminate this CT need for MRgFUS and MR/PET.
The preliminary implementation of preconditioned TFI can be further developed with the following considerations. First, the current choice of the preconditioner P is empirical and performance may be improved by selecting a patient specific value or model. Further work might focus on determining P from physical or statistical models of the susceptibility distribution for better approximation of the covariance matrix Γ. A prior QSM estimation, which can be calculated very fast by methods such as SSQSM, might also be helpful in modeling this distribution and constructing the preconditioner. Second, the construction of our preconditioner requires R2* map to extract intracerebral region with strong susceptibility, such as ICH. If the R2* information is not available, especially when only single-echo image is acquired, the preconditioned TFI might be less effective in suppressing ICH-based shadow artifact. In future work, we will focus on incorporating other contrast (magnitude, T1-weighted or T2-weighted image) into constructing the preconditioner. Third, for whole head QSM, the skull has very short T2*, such that the multi-echo GRE sequence fails to acquire skull MR signal for field estimation. In future work, the use of an ultra-short echo time (UTE) pulse sequence to provide phase information in the skull will likely lead to a more accurate susceptibility map. The magnitude signal from UTE sequence might also be used for differentiating skull and air-filled sinus, hence enabling more effective preconditioning by assigning different weight PB for these regions. Finally, the proposed TFI method uses a linear signal model but can be extended to a non-linear formulation similar to that of, in order to bypass the phase unwrapping step and improve the noise modeling. However, caution is needed because the non-linear and non-convex objective function has multiple local minimum. Given the fact that the desired solution is fairly large outside the ROI, such as that of air, the all-zero initial guess used in this work may cause the solver to converge to an incorrect local minimum. Future work will be focused on construction of an initial guess robust against this problem. The current linear solution may provide such an initial guess.
In summary, we describe a preconditioned total field inversion (TFI) algorithm for QSM reconstruction. This eliminates the need for sequential background field removal and local field inversion, and it achieves a higher accuracy at the boundary of the susceptibility map. We also introduce a preconditioner to improve convergence. We show that the proposed preconditioned TFI suppresses streaking artifacts in QSM of intracerebral hemorrhage. TFI is also capable of mapping the susceptibility of the entire head including brain tissue, skull, air-filled sinus and subcutaneous fat.
In this section, we present an improvement on constructing the structure prior by taking into account of the structural edge orientation.
Using a mathematical analysis of the field-to-susceptibility inverse problem, it has been shown that streaking artifacts in QSM are originated from the dipole-incompatible part of the field. This incompatible part includes noise, discretization induced error, anisotropic sources or chemical shift. Furthermore, it is shown that streaking-free solutions can be obtained from the dipole-compatible part of the field only. However, noise is omnipresent in practice, discretization is indispensable for numerical computation, and anisotropic sources cannot be measured reliably in a single-angle acquisition. Therefore, regularization in image space plays a crucial role in computing solutions with minimal streaking.
Among others, morphology enabled dipole inversion (MEDI) is an approach capable of enhancing tissue structure by imposing spatial consistency (edge information) between a candidate susceptibility map and the magnitude image associated with the susceptibility map. To be specific, the structural prior is incorporated into total variation (TV) regularization as an isotropic matrix. This aims at 1) preserving common edges between the magnitude image and a candidate susceptibility map—tissue structure; and 2) penalizing strong gradients elsewhere—streaking artifacts. Therefore, its performance largely depends on how to design such a structural prior (matrix). This paper focuses on the design of this structural prior that takes into account orientation information which was previously missing in isotropic weighting.
The remainder of this paper is organized as follows: In THEORY, we consider anisotropic weighting (AW) for a structural prior that takes into account orientation consistency. In METHODS and RESULTS, we demonstrate the effectiveness of AW in two in vivo brain MRI datasets where COSMOS or the component of the susceptibility tensor X33 are available as reference. In DISCUSSION and CONCLUSION, we conclude with some discussion on our results. MATLAB source codes will be available at http://weill.cornell.edu/mri/.
In a spatially continuous setting, the problem of regularized dipole inversion within a region of interest Ω is written as the following variational problem:
where the susceptibility distribution χ:(Ω⊂3)→
(L1-integrable) is estimated from the measured local field b:Ω→
. The admissible solution space BV(Ω;
) of bounded variation (BV) allows the existence of jumps (sharp edges). The map w: L1→L1 in the data fidelity term compensates for the non-uniform phase noise (SNR weighting). Here, d and* are the unit dipole kernel and the convolution operation, respectively. The regularization term Reg(χ) is designed in such a way that it effectively penalizes streaking artifacts originated from the dipole-incompatible part of the field b. In MEDI, for instance, the regularization term is formulated as follows:
where δx(r), δy(r), and δz(r) are edge indicators for the x-, y-, and z-directions at r=(x,y,z), respectively. That is, δx(x0, y0, z0)=1 if there is an edge along the x-axis at r0=(x0, y0, z0); otherwise δx(r0)=0. Likewise, δy(r) and δz(r) are defined. Since these edge indicators are computed from the magnitude image associated with the unknown susceptibility map under the assumption that edges in the two images coexist, the regularizer Reg(χ) penalizes only the regions where tissue structure (edges) is not expected. Note that λ>0 is a tuning parameter that balances data fidelity and regularization.
In order to take into account edge orientation information, an anisotropic weighting can be defined as follows:
where ξ(r) is the edge indicator δ:Ω→, r
δx(r)∨δy(r)∨δz(r) multiplied by the gradient vector of the magnitude image Mag(r) at r as follows:
Note that ∨ denotes the Boolean OR operation. The matrix P⊥ξ is the rank-2 orthogonal projector that eliminates the component in the direction of ξ. In other words, we modify the regularization term such that it penalizes the components of the susceptibility distribution gradient that are orthogonal to the magnitude image gradient, as illustrated in
Not only is the data fidelity term convex, but also both regularization terms Reg(χ) and Reg⊥(χ) are convex in χ. Therefore, one can use some well-established algorithms for this convex optimization problem such as the Gauss-Newton conjugate gradient method, the split-Bregman method, and the primal-dual algorithm to efficiently minimize the energy [12]. In this paper, we use the Gauss-Newton conjugate gradient method with the Sobolev space assumption, which solves the Euler-Lagrange equations associated with the first variation of the energy functional in [12]. The first variation of Reg(χ) is derived as follows:
where the gradient operator ∇ should be understood in the distributional sense in BV. Note that the above first variation [16] is not well-defined in BV but is in the Sobolev space, which is the assumption behind the Gauss-Newton conjugate gradient method.
The first variation of Reg⊥(χ) is derived as follows:
Likewise, this first variation is defined in the Sobolev space.
Here we describe experimental details to evaluate the implementation of anisotropic weighting in TV regularization with respect to isotropic weighting.
Two in vivo brain MRI data in 2 healthy subjects were used for experimental validation of the effects of anisotropic weighting using COSMOS or χ33 as reference. These data sets were obtained from the Cornell MRI research group (http://weill.cornell.edu/mri/pages/gsm.html) and from the QSM 2016 reconstruction challenge (http://www.neuroimaging.at/qsm2016/pages/qsm-challenge.php). For imaging parameters and phase preprocessing, we refer the reader to the websites and references therein.
Numerical conditioning was set to ∈=10−6; the maximum number of inner iterations (nested CG loop) was set to 100; the maximum number of outer iterations was set to 100; the tolerance for the nested CG loop was set to 0.01; and the relative update norm ∥χk+1−χk∥/∥χk∥, where χk is a solution at the kth outer iteration, was set to 0.01. Note that the stopping criteria are made tighter than default settings to make sure convergence—typically the maximum number of outer iterations is 10 and the relative update norm is 0.1. The algorithm was run in Matlab R2016a on an Intel Core i7 3.30 GHz CPU with 64 GB of memory.
In isotropic weighting, the edge indicators δx, δy, and δz are determined as the following procedure: 1) Compute ζx(r)=|∂xMag(r)|, ∂y(r)=|∂yMag(r)|, and ∂z(r)=|∂zMag(r)| for r∈Ω; 2) define the edge indicators as characteristic functions of a parameter t as follows: δx(r)={ζ
{ζ
{ζ
Then, the resulting binary indicators δx(r), δy(r), and δz(r) are used for [13] and [16].
In anisotropic weighting, we need to determine ξ(r)=(ξ1(r),ξ2(r), ξ3(r))T. This boils down to computing δ:Ω→ as defined in [15]. Once δx, δy, and δz are computed by the procedure described above in isotropic weighting, we determine δ as δx(r)∨δy(r)∨δz(r) for r∈Ω. These two weighting methods were visualized in
As the values of λ/2, 1000 and 235 were empirically chosen using the discrepancy principle for the Cornell MRI research lab data and the QSM 2016 reconstruction challenge data, respectively. Then, in the range of 700 to 1350 (700, 750, . . . , 1350) for the Cornell MRI research group data and in the range of 164.5 to 305.5 (164.5, 176.25, . . . , 305.5) for the QSM 2016 reconstruction challenge data, the MEDI algorithm with anisotropic and isotropic weighting was performed.
The following evaluation metrics were used for QSM accuracy assessment: RMSE (root mean square error), HFEN (high frequency error norm), and SSIM (structure similarity index). We used MATLAB codes for RMSE and HFEN available at http://www.neuroimaging.at/qsm2016/pages/qsm-challenge.php, and for SSIM we used the MATLAB function (SSIM) in Image Processing Toolbox with default settings. In addition, we used what we call the degree of alignment to measure how well aligned two images are, which is defined as follows: Let u and v be scalar-valued images (e.g., functions in BV(Ω; )); and let nu, nv denote vector fields of surface normals to u and v, respectively. Suppose that θn
where |⋅| denotes the Lebesgue measure in 3, i.e., volume. Note that the quantity DoA(u, v) is bounded below by 0 (when the two vector fields are perpendicular almost everywhere) and is bounded above by 1 (when the two vector fields are parallel—either 0 rad or π rad). Therefore, the quantity DoA(u, v) measures how well aligned the two images u and v are in the range of [0,1] (higher indicates better alignment).
To assess QSM accuracy based on mean voxel values in various brain regions, 19 circular regions of interest, each measuring 14.5 mm3, were drawn on selected brain slices by a neuroradiologist with more than 20 years of experience, in the striatum (caudate nucleus, putamen, globus pallidus), thalamus, cortical gray matter, splenium of corpus callosum, frontal white matter, occipital white matter, and cerebrospinal fluid (CSF) in lateral ventricles. The mean voxel value for each ROI was calculated in MATLAB and recorded relative to the CSF mean value. For various values of regularization parameter, λ, the slope and coefficient of determination, R2, were recorded for the linear regression between the MEDI derived ROI values and the COSMOS or χ33 derived ROI values. These fitting parameters were compared between isotropic or anisotropic weighting for MEDI.
For the Cornell MRI research group data and for the values of λ/2 (700 to 1350), the CPU running times were approximately 38 min and 44 min for anisotropic and isotropic weighting, respectively. For the QSM 2016 reconstruction challenge data and for the values of λ/2 (164.5 to 305.5), the CPU running times were approximately 10 min for both anisotropic and isotropic weighting.
RMSE (lower is better), HFEN (lower is better), SSIM (higher is better), and the DoA (higher is better) were shown in
QSM accuracy was assessed in different methods (isotropic and anisotropic) via linear regression based on voxel values within ROIs between the methods and reference standards (COSMOS or χ33).
We have considered anisotropic weighting as a new structural prior in the QSM inverse problem. As opposed to isotropic weighting, anisotropic weighting imposes orientation consistency between a candidate susceptibility map and the magnitude image that is associated with the susceptibility map. In experiments, we have shown improvements in QSM accuracy in terms of both quantitative measures and image quality within a wide range of regularization parameter.
The notion of anisotropic weighting considered in this paper was widely studied in a seminal monograph on anisotropic diffusion in image processing and in the context of vector-valued image processing whose main focus is to retrieve common edges across the channels (e.g., RGB) and suppress random noise. Recently, a similar work has been studied in image processing in the name of directional total variation and another work has been published in the context of multicontrast MRI for fast imaging where structural information in T1 weighted images were used for the reconstruction of T2 weighted images with sparse samples. The proposed method is in line with these works; but is new in quantitative susceptibility mapping.
In
Therefore, when θ(r0)=π/4 and 5π/4, there is no penalty (maximally allowed); when θ(r0)=3π/4 and 7π/4, penalty is maximized (maximally disallowed). In isotropic weighting, on the other hand, the following is regularized:
where δx(r0) and δy(r0) are edge indicators for x and y directions, respectively. Therefore, we have |sin θ(r0)|, |cos θ(r0)|, and 0 penalties for the x, y, and both x and y directional edges, respectively. Clearly, this does not take into account orientation consistency—for instance if both x and y directional edges are detected, every possible direction of ∇χ(r0), i.e., θ(r0) is allowed with no penalty. This additional “adaptive penalization” in anisotropic weighting appears to have suppressed such visual artifacts originated from orientation inconsistency.
The CPU time gaps between anisotropic and isotropic weighting for the Cornell MRI research group data and the QSM 2016 reconstruction challenge data were approximately 6 min and few seconds, respectively. Anisotropic weighting was faster on the Cornell MRI research group data and was as fast as isotropic weighting on the QSM 2016 reconstruction challenge data. Anisotropic weighting involves 3×3 full matrix multiplications as opposed to 3×3 diagonal matrix multiplications for isotropic weighting leading to few more seconds at each outer iteration of the Gauss-Newton conjugate gradient method. However, oftentimes anisotropic weighting converged faster with less number of outer iterations, which explains the fact that anisotropic weighting is as fast as isotropic weighting.
To devise a structural prior that takes into account the biophysical model behind QSM is the cornerstone of the success for regularization-based inversion methods. As seen throughout the paper, the anisotropic weighting approach outperforms isotropic weighting and its implementation is very simple. We believe that the proposed approach will be an important step towards an optimal field decomposition method; albeit unknown yet.
In summary, incorporating anisotropic weighting into the regularization term of MEDI improves QSM accuracy by taking into account orientation consistency which is missing in isotropic weighting.
In this section, we describe the primal dual solver for the QSM problem, which would improve both accuracy and speed for evaluating the derivative of the L1 norm term.
Bayesian MAP estimation has been used increasingly in MRI to reconstruct images from incomplete noisy data with a likelihood to model noise and a prior probability to compensate for undersampling. Its application in image reconstruction often leads to minimization problems that can be commonly solved by using the Gauss-Newton conjugate gradient (GNCG) algorithm. This is the case for quantitative susceptibility mapping (QSM), where a structural prior is used to identify a solution with minimal streaking (morphology enabled dipole inversion or MEDI). The structural prior imposes the spatial consistency of edges (gradients) between a known tissue image and the targeted tissue susceptibility map.
The goal of this paper is to investigate the effects of the specific implementation of this gradient L1-norm prior term on QSM quality. First, typical GNCG implementations employ numerical conditioning that slightly modify this prior term smoothing the L1-norm, thus solving an approximated energy function. While the effects of this conditioning have been investigated in the context of image denoising, an analysis for QSM has not been performed yet. An alternative method, the primal-dual algorithm, will be used for comparison. The latter algorithm has been used to compute generalized total variation in Bayesian MRI and in QSM. Second, a central difference scheme has been used to implement the gradient operator in the prior term. This can lead to checkerboard artifacts. An alternative method, a forward difference scheme, will be used for comparison. Accordingly, this paper will evaluate the primal-dual formulation of QSM using the Chambolle-Pock algorithm with a forward gradient implementation and compare this with the GNCG method with a central gradient implementation.
Here, we outline a spatially continuous formulation for the QSM energy minimization problem to describe the dual definition of the anisotropic total variation, instead of the approximate derivative used in GNCG. Details on the primal-dual algorithm can be found in the MEDI toolbox (http://weill.cornell.edu/mri/pages/qsm.html).
Currently, QSM has been formulated as a discrete energy minimization problem to estimate the underlying susceptibility distribution. Its continuum limit can be expressed by the following variational model:
Here, the susceptibility function u: Ω→ (L1-integrable) is estimated from the measured local field b: (Ω⊂
3)→
(Ω is a Lipschitz domain). The solution space BV(Ω;
) of bounded variation (BV) allows the existence of jumps (sharp edges). The objective consists of two energy terms, data fidelity and regularization, as in the discrete formulation. The data fidelity term involves a linear operator T defined by Tu(r)=(d*u)(r), r∈Ω, where d and* are the unit dipole kernel and the convolution operation, respectively. Let L1 be the space of absolutely integrable (L1-integrable) functions. The map W: L1→L1 compensates for the non-uniform phase noise (SNR weighting). The regularization term imposes spatial regularity using anisotropic weighted total variation (TV) with an edge map M: [L1(Ω;
)]3→[L1(Ω;
)]3 that is extracted from the corresponding magnitude image. The gradient operator D is understood in the sense of distribution. The tuning parameter λ>0 balances data fidelity and regularization, which can be determined by using the discrepancy principle or the L-curve method, etc. Based on this spatially continuous formulation of MEDI, the existence of a minimizer can be established using the BV deblurring theory. The formulation leads us to functional minimization where the susceptibility u is a function of r. Then, the susceptibility denoted by χ commonly used in QSM literature is given by discretizing u. Our interest here is to compute such a minimizer, and we start with the GNCG algorithm.
The GNCG method assumes that u is differentiable (technically u is assumed to be in the Sobolev space W1,1), which allows us to derive the Euler-Lagrange equation associated with the energy functional as follows:
where (WT)* is the adjoint of WT. The Euler-Lagrange equation, however, exhibits degeneracy when |(r)
u(r)|=0 so that numerical conditioning is required for computational purposes: |
(r)
u(r)| is replaced with
and 3 and a small ∈>0. As a consequence, this conditioning deviates from the definition of TV for L1-integrable functions where image edges are modeled as the set of discontinuities. Therefore, the algorithm solves a slightly different objective that involves the conditioning parameter ∈.
Similar to the distributional derivative for a function with discontinuities, we can derive a generalized definition of anisotropic weighted TV in the space L1 of integrable functions
∫Ω|MDu|1=sup{−∫Ωu div(M*ξ)dr|ξ∈Cc1(Ω,3),∥ξ∥∞≤1,∀r∈Ω} [20]
where Cc1 is the set of continuously differentiable functions with compact support, ∥ξ(r)∥∞=max{|ξx(r)|, |ξy(r)|, |ξz(z(r)|} and M* is the adjoint of M. ξ is the so-called dual variable with the constraint of pointwise maximum norm, the divergence of a vector is denoted by div, and the supremum of a set is denoted by sup. Note that the discontinuities are allowed in u in the above dual formulation. Since
is self-conjugate and (αH)*=αH*(⋅/α) for H:L2→(−∞, +∞] and α>0, the data term can be dualized as follows:
where ⋅, ⋅
L
). Therefore, the minimization problem [18] can be rewritten using [20] and [21] resulting in the so-called primal-dual formulation:
where the indicator function tP(ξ) of a set P⊂Cc1 is defined as
with the convex set
P={ξ∈C
c
1(Ω;3):∥ξ(r)∥∞≤1,∀r∈Ω}. [23]
This is a saddle-point problem with a convex constraint. Under the assumption that there exists a saddle-point satisfying solution [22], a provably convergent algorithm can be applied to compute such an optimal point. Among others, we make use of the fast primal-dual (PD) algorithm proposed by Chambolle and Pock. In the following, we have provided a discrete setting for [22] and implementation details on the PD algorithm to compute an optimal susceptibility map. Note that GNCG now can be viewed as the primal approach.
To compute an optimal solution of the saddle-point problem [22], first we provide a discrete setting. Let us consider a three-dimensional regular Cartesian grid of size Nx×Ny×Nz as follows:
{(ihx,jhy,khz):1≤i≤Nx,1≤i≤Ny,1≤i≤Nz},
where hx, hy, and hz denote the voxel size, and (i, j, k) the voxel location. Define a scalar product on H=N
N
N
For u∈H, the discrete gradient ∇u is defined by
The gradient vector ∇u is a vector in the finite dimensional real vector space K=H×H×H where a scalar product is given as follows:
Notice that the forward difference scheme with the Neumann boundary condition is used for the discrete gradient. Now, we define a linear map M:K→K such that
ξi,j,k(Mξ)i,j,k=(Mi,j,kxξi,j,kx,Mi,j,kyξi,j,ky,Mi,j,kzξi,j,kz)
for all i, j, and k. Then, the following discrete dual definition of [20] is straightforward
∥M∇u∥1=sup{u,∇*M*ξ
H:ξ∈K,∥ξi,j,k∥∞≤1,∀i,j,k}
where ∇*:K→H and M*:K→K are the adjoints of ∇ and M, respectively; i.e., M∇u, ξ
K=
u, ∇*M*ξ
H. The L∞-norm ∥⋅∥∞ is defined as
∥ξi,j,k∥∞=max{|ξi,j,kx|,|ξi,j,ky|,|ξi,j,kz|}.
From the relationship u, ∇*M*ξ
H=
M∇u, ξ
K and the discrete setting for M∇u, the following discretization of the adjoint operator ∇*ξ is easy to derive:
Note that the adjoint operator of the discrete gradient is the negative divergence, i.e., ∇*M*ξ=−div(M*ξ), which is discretized by the backward difference scheme with the Dirichlet boundary condition.
The function spaces Cc1 and L2 are replaced by K and H, respectively; and, the convex set [23] is discretized as follows:
P={ξ
i,j,k
∈K:∥ξ
i,j,k∥∞≤1,∀i,j,k}.
Then, a discrete formulation of [22] is now straightforward. This discrete primal-dual formulation can be efficiently solved by a provably convergent saddle-point solver. Here, we apply a first-order primal-dual algorithm whose convergence was proved by Chambolle and Pock with an overrelaxation step. The algorithm is basically read as alternating gradient descent with respect to the primal variable u followed by the proximal operator associated with W(Tu−b), v
H and gradient ascent with respect to the dual variables ξ and v followed by the proximal operators associated with ιP(ξ) and
Let us define G as W(Tu−b), v
H, then the proximal operator of the scaled function τG, τ>0 is given as
u=proxτG(ũ)⇔u=ũ−τT*W*v⇔u=ũ−τ−1(
⊙
(W*v)),
where (⋅) and
−1(⋅) denote the fast Fourier transform (FFT) and inverse FFT, respectively.
is the complex conjugate of
(d) and the operation ⊙ is the pointwise multiplication. Note that at each iteration, only one FFT and one inverse FFT are required to perform for this operator;
is computed in advance. On the other hand, the proximal operators of the scaled function σιP(ξ) and
are given as
for all i, j, and k. Finally, we describe the primal-dual algorithm as follows:
for n=0, . . . , N−1. Again, to update v at each iteration one FFT and one inverse FFT are required. As initialization, we choose (u0,ξ0)∈H×K and set ū0=u0. The primal and dual step sizes τ and σ are chosen such that τσ≤1/∥M∇∥2 to ensure the convergence of the algorithm. Since 0≤Mi,j,kx≤1, 0≤Mi,j,ky≤1, and 0≤Mi,j,kz≤1, one has:
Here we describe experimental details to implement and evaluate the primal-dual QSM with the forward gradient.
A simulated brain, a gadolinium (Gd) phantom, in vivo brain MRI data in 4 healthy subjects, and ex vivo brain blocks were used for experimental validation of the implementations using as reference the ground truth when available (simulation and Gd concentration) and COSMOS for the healthy subjects and brain blocks. The simulated brain and one set of the in vivo brain MRI data were downloaded from http://weill.cornell.edu/mri/pages/gsm.html. The rest of the in vivo data were extracted from. For phase preprocessing, a temporal unwrapping of the phase was performed in each voxel followed by a weighted least-squares fit of the temporally unwrapped phases in each voxel over TE. Then, a magnitude map guided spatial unwrapping algorithm was applied. The local field is then extracted by the projection onto dipole fields (PDF) algorithm. Note that the white matter of the simulated brain is piecewise smooth, and the rest of the brain (gray matter, etc.) is piecewise constant.
Excised brain blocks from three postmortem control brains and three postmortem Parkinsonian brains (Yale-New Haven hospital, New Haven, Conn.) were fixed and embedded in 1% agarose gel. The brain blocks were scanned using a 3T clinical MRI system (GE Healthcare, Milwaukee, Wis.) using a standard 8 channel head coil. Multiple-echo 3D gradient echo (GRE) data were acquired using the following parameters: 0.4 mm isotropic resolution, TE=4.3 ms, ΔTE=4.8 ms, #TE=11, TR=74.2, BW=62.5 kHz, FA=20°. Three in-plane orientations (Δϕ=60°) were acquired. For phase processing, we used a temporal fitting method and a spatial phase unwrapping algorithm.
Let hx, hy, and hz denote the discretization widths, or voxel size, and (i,j,k) denote the indices of the discrete locations (ihx, jhy, khz) in the (volumetric) image domain; see the above for implementation details. Then, a discrete version of the anisotropic TV of u based on the central difference scheme for the gradient operator is defined by
and another discrete version of the anisotropic TV of u based on the forward difference scheme for the gradient operator is defined by
Then, a discrete version of the anisotropic weighted TV is straightforward. GNCG is implemented using the central difference scheme (GNCG-C) and the forward difference scheme (GNCG-F). PD is implemented using the forward difference scheme.
For GNCG, numerical conditioning was set to ∈=10−6, the maximum number of inner iterations (nested CG loop) was 100, the maximum number of outer iterations was 20, and the tolerance for the nested CG loop and for the relative update norm ∥uk+1−uk∥/∥uk∥, with uk the solution at the kth outer iteration, were both set to 1%. For the PD algorithm, the relative update norm was set to 0.0001% and we set τ=1/(σ∥M∇∥2) to ensure convergence, where σ=∥M∇∥2 for ex vivo brain slices and σ=∥M∇∥4 for the rest of the dataset. These algorithms were implemented in Matlab R2015b on an Intel Core i7 3.30 GHz CPU with 64 GB of memory. Additionally, PD was implemented in MATLAB (unoptimized) that supports GPU computing and was run on NVIDIA GPUs (NVIDIA GeForce GTX 750 Ti and TITAN X).
GNCG-C, GNCG-F, and PD were run with various parameter settings (∈ for GNCG and τ and σ for PD) on a simulated brain. Stopping criteria for these algorithms were set to be tighter than previously described to allow for investigation of the convergence behavior of each algorithm. To this end, we computed the relative norm defined by ∥uk−u*∥2/∥u*∥2, where u* is the ground truth, at each iteration k during optimization.
The COSMOS image for each in vivo and ex vivo data set was loaded into the ITK-SNAP image processing program and regions of interest (ROIs) were drawn. ROIs for the in vivo data included the cortex, deep gray matter structures (red nucleus, substantia nigra, caudate nucleus, putamen, globus pallidus) and white matter regions (splenium of corpus callosum, frontal white matter, occipital white matter). ROIs for ex vivo data included the cortex, putamen, globus pallidus, red nucleus, substantia nigra and subthalamic nucleus.
Linear regression was performed by fitting the regression line, X=k1XCOSMOS+b, where X was the set of ROI mean values for all subjects for either GNCG-F or PD and XCOSMOS was the set of corresponding COSMOS values. The slope k1, intercept b, and coefficient of determination R2 were recorded for each regression. To further assess agreement between COSMOS and the two reconstruction methods, Bland-Altman plots were drawn to determine the mean difference (bias) and 95% limits of agreements.
Our results are organized according to convergence analysis, gadolinium phantom, ex vivo brain slabs, and in vivo imaging.
The number of flops required for GNCG and PD algorithms per iteration is almost the same. Here, iteration refers to one cycle of the nested CG loop for GNCG (#operations=2 FFTs+2 iFFTs+3 dot products+2 SAXPYs+2 GAXPYs) and one cycle of the algorithm for PD (#operations=2 FFTs+2 iFFTs+4 SAXPYs+2 GAXPYs). Note that for GNCG additional operations (1 SAXPY and 1 GAXPY) are added from the outer loop once the nested CG loop is terminated. For GNCG, the maximum number of inner iterations per outer loop was 100, with actual number of inner iterations varying from outer loop to outer loop and from dataset to dataset (50˜100).
The CPU running times of each algorithm per iteration were approximately 0.79 sec, 0.80 sec, 0.75 sec, 0.75 sec, 0.75 sec, and 0.75 sec for GNCG-C (ε=10−6), GNCG-C (ε=10−8), GNCG-F (ε=10−6), GNCG-F (ε=10−8), PD (L3), and PD ((L3+L4)/2), respectively.
The accuracy of GNCG-F was: for λ/2=150, linear regression slope=1.06, R2=0.79; for λ/2=300, slope=1.06, R2=0.92. The accuracy of PD was: for λ/2=150, slope=1.04, R2=0.87; for λ/2=300, slope=1.04, R2=0.81. The CPU running times (size of the Gd phantom=130×130×116) of each algorithm per iteration were approximately 0.36 sec for both GNCG-F and PD, and the GPU (GTX 750 Ti) running times of PD per iteration were approximately 0.03 sec. The total number of iterations for GNCG-F was 958 (with 10 outer iterations). The total number of iterations for PD was 800.
In this paper, we show the influence on the convergence rate and accuracy of various implementations of the regularization term in QSM. The numerical conditioning of this term necessary for the use of the GNCG solver leads to slower convergence as compared to the continuous formulation of QSM as a saddle-point problem using the primal-dual (PD) algorithm. The use of a central difference approximation for the gradient in the regularization term can lead to checkerboard artifacts that are greatly reduced when using a forward difference approximation.
The Chambolle-Pock (primal-dual) algorithm has been applied to total generalized variation (TGV) in the context of QSM, and the split-Bregman method was alternatively used for fast QSM. Our contributions here include the following aspects: First, unlike simplified data fidelity terms (least squares) used in, we consider the MEDI problem where SNR weighting is involved (weighted least squares), which is necessary to reduce noise error in QSM but significantly increases computational cost because of the addition of two 3D Fourier transforms per CG iteration and PD iteration. Second, we employ the weighted anisotropic TV, rather than TGV, to enforce structural consistency between known tissue structure and targeted susceptibility map (morphology enable dipole inversion, MEDI). This physical MEDI has been quantitatively compared with other existing methods on clinical data to be advantageous over generic sparsity defined by unweighted TV or TGV. Third, we provide the perspective of variational methods for QSM.
Note that the QSM problem in [18] is different from the image deblurring problem considered in in the following aspects: First, as opposed to a conventional blur operator in image deblurring, the forward operator in QSM, namely the dipole kernel (with SNR weighting), does not decay to zero away from the origin in k-space but vanishes on the conic surface spanned by the magic angle (≈54.7°) from the main magnetic field. This fundamental difference in the forward operators leads to somewhat peculiar spectral properties in QSM, which was investigated in terms of the discrete Picard condition. Second, the role of TV regularization is different. Streaking artifacts can be avoided by extracting the dipole-compatible part of the field only. Although such a decomposition method is not known, the anisotropic weighted TV in [20] turns out to be able to effectively penalize streaking artifacts while preserving the underlying tissue structure. The TV regularization term in image deblurring, on the other hand, plays a role in edge-preserving smoothing while reducing additive white Gaussian noise.
Our work here illustrates the conditioning dependent convergence and accuracy in the standard GNCG, which are inferior to PD (
The faster convergence of PD over GNCG (
The suppression of checkerboard artifacts by the forward difference approximation for the gradient operation can be seen in the difference between GNCG-C and GNCG-F in
and the forward difference approximation defined by
For the nested 3×3 square of the matrix, we have Jc=0 and Jƒ=12. Indeed, there is no penalization by the central difference scheme, whereas there clearly exist constant fluctuations in u along the both x and y axes, which are captured by the forward difference scheme.
The error in GNCG may be understood by Huber norm analysis. Under the Sobolev assumption, the notion of numerical conditioning in the Euler-Lagrange equation of MEDI has a close relationship with the Huber-norm defined by
where α>0 is a parameter that determines a trade-off between quadratic regularization for small ranges of α and anisotropic TV regularization for larger ranges of α. Replacing the integrand in the regularization term in [18] with the Huber-norm, we can derive its Euler-Lagrange equation when |MDu|1≤α as follows:
A time-marching method can be applied to solve this equation; the overall task involves solving a diffusion equation whose diffusivity constant is given as ½α. Indeed, such diffusion processes take place in the regions where |MDu|1≤α. Now, consider the Euler-Lagrange equation [19] with the conditioning parameter ∈ in its denominator, i.e., as described in the THEORY section. Applying a time-marching method to solve this equation, we notice that if |
(r)
(r)|<<1, then
≈2α for some ∈ and α. Hence, the smooth regions shown in
QSM accuracy may be further improved by recent technical advances in nonsmooth convex optimization. Energy minimization methods based on sparsity-inducing norms such as TV-seminorm have been a viable solution in a broad range of ill-posed imaging problems. It is, however, preferable to work with an energy function that faithfully captures the underlying physical nature of the imaging data. Otherwise, computing an optimal solution for a modified energy function may lead to an undesirable solution, as seen, for example, in the role of proper noise modeling for the data fidelity or prior term. The discrepancies between MEDI and COSMOS observed in ex vivo and in vivo data analyses may also be explained by model error in ignoring white matter susceptibility anisotropy and by mis-registration error.
In summary, we have shown the influence of specific implementation of the regularization term in QSM as it relates to convergence rate and accuracy. The slower convergence rate and accuracy associated with the numerical conditioning required for standard GNCG is substantially improved using a primal-dual algorithm that no longer requires the modification of the energy function in QSM. In addition, we have shown that the forward difference scheme for implementing the gradient operator suppresses the checkerboard artifacts caused by the central difference scheme. Over various types of data from phantom to in vivo human brain, the primal-dual formulation with forward gradient has shown improved convergence and accuracy over the current GNCG.
In this section, we describe a QSM improvement that provides a uniform cerebrospinal fluid (CSF) in ventricles for automatic zero reference, as well as general shadow artifacts reduction.
Current quantitative susceptibility mapping (QSM) methods quantify magnetic susceptibility with respect to a certain level, known as the zero susceptibility reference. For brain parenchymal QSM, CSF is widely used as the zero reference because of its chemical resemblance to pure water. However, in current QSM, there is often a non-uniformity in the CSF, which increases the uncertainty of the reference susceptibility. Literature suggests that CSF flow and white matter anisotropy might contribute to this variation in CSF susceptibility.
Considering that streaking artifacts belonging to the null space of the field-to-source inverse problem can be penalized using a tissue structure prior in current Bayesian QSM (morphology enabled dipole inversion, MEDI), we regard CSF susceptibility inhomogeneity also as belonging to the null space, and we propose to use a prior that enforces the susceptibility homogeneity inside the CSF. In this paper, we present preliminary data demonstrating that minimal susceptibility variation in CSF can be achieved for a consistent and automated zero reference of QSM, which is referred to as QSM0. Furthermore, QSM0 brings about additional image quality improvement.
The cost function in the optimization problem for QSM in previous literature (MEDI) consists of a data fidelity term corresponding to the Gaussian noise in the complex MRI data, and a regularization term reflecting the structure prior information regarding the magnetic susceptibility distribution, which is the sparse discrepancy between known tissue edges and edges in the susceptibility map:
with χ the susceptibility map, * the convolution with the dipole kernel d, w the noise weighting, b the magnetic field experienced by tissue with Lorentz correction or the local field or simply the magnetic field, ∇ the gradient operator and MG the binary edge mask derived from the magnitude image. The proposed QSM0 introduces an additional regularization term into the cost function to enforce a known uniform susceptibility in a region of interested tissue, solving the following optimization problem:
Here the additional term to MEDI (Eq. 24) is an L2-regularization, where CSF is the interested tissue with a known uniform susceptibility, and MCSF is the region of interest (ROI) of lateral ventricular CSF with a mean susceptibility
We propose an automated procedure to determine the lateral ventricle mask MCSF using the brain ROI mask M and R2* map. CSF with its chemical composition mostly made of water is known to have a characteristic feature of very low R1, R2 and R2* values. Accordingly, CSF can be segmented out from anatomic images, such as the magnitude images from gradient echo MRI used in QSM or R2* map generated these magnitude images, based on a characteristic feature of CSF on anatomic images, such as uniform low R2* value, which allows thresholding. Multiple regions are collected by applying connectivity. A specific implementation includes the following steps:
while CSF presumably has a smaller R2* than the threshold R.
where N is number of voxel in M.
divide Mc∩MR{Mc1∪Mc2∪Mc3}.
divide
into connected components Mi (6-neighbour) and merge all components overlapping with McCSF: MCSF{∪Mi|Mi∩McCSF≠Ø}.
The proposed algorithm introduces two new parameters: λ2 and R, which will be determined in a numerical simulation discussed later.
The problem (Eq. 25) is solved using Gauss Newton Conjugate Gradient. As a final step,
We constructed a numerical brain phantom to simulate the CSF inhomogeneity related to anisotropic susceptibilities. Matrix size was 256×256×126, and voxel size was 1×1×1.5 mm. The susceptibility tensor was originally constructed from 12-orientation human brain scans using susceptibility tensor imaging (STI). Lateral ventricles were manually drawn (Y.Y., 4 years of experience), and regarded as the region of CSF. Under the assumption that CSF is homogeneous hence not a source for susceptibility anisotropy, susceptibility was set to 0 for CSF in each element of the susceptibility tensor. Then the synthetic local tissue field was generated with B0=[0, 0, 1] in two ways: first, tissue field ƒaniso was generated with χ13, χ23, χ33 using the anisotropic forward model:
where F and X are the Fourier transform of f and χ respectively; Second, tissue field ƒiso was generated solely with χ33 using the isotropic forward model,
Finally, Gaussian noise (SNR=100) was added to both fields. MEDI and QSM0 were applied on both fields. The root mean square error (RMSE) between the reconstructed QSM and χ33, Err=∥χ−χ33∥2, was calculated for each method. λ1 was determined by minimizing Err, and used for both MEDI (Eq. 24) and QSM0 (Eq. 25). λ2 and R in QSM0 were also determined from minimizing Err.
The brains of 8 patients with Multiple Sclerosis (MS) were imaged using multi-echo GRE at 3T (GE, Milwaukee, Wis.). All studies in this work were approved by our Institutional Review Board. Imaging sequences were: T2w fast spin echo (Flip angle=90, FOV=24 cm, TE=86 ms, TR=5250 ms, slice thickness=3 mm, acquisition matrix=416×256×68) and 3D T2*w spoiled multi-echo GRE (Flip angle=20, FOV=24 cm, TE1=4.3 ms, TR=57 ms, #TE=11, ΔTE=4.8 ms, acquisition matrix=512×512×68, voxel size=0.47×0.47×2 mm3, BW=±62.5 kHz, flow compensation with monopolar readout, parallel imaging factor R=2 and a total scan time of ˜10 min). Nonlinear field map estimation followed by graph-cut-based phase unwrapping was applied to estimate the total field, then projection onto dipole fields (PDF) was used for background field removal in order to obtain the local field. Next, MEDI and QSM0 were performed on the same local field. The standard deviation of susceptibility in ventricular CSF was measured for both methods. For lesion susceptibility measurement, one neuroradiologist (Y.Z., 4 years of experience) drew the ROIs of each lesion on the T2w images, which were already registered onto magnitude images of the T2*w scan. Reference regions for lesion susceptibility were determined by drawing an ROI in the normal appearing white matter (NAWM) on the contralateral side of the identified lesion. The difference of mean susceptibilities between the lesion and NAWM were measured and analyzed using linear correlation and Bland-Altman plots to quantify the agreement between MEDI and QSM0. Analyses were also performed on ROIs for the in vivo data included the cortex, deep gray matter structures (red nucleus, substantia nigra, caudate nucleus, putamen, globus pallidus) and white matter regions (splenium of corpus callosum, frontal white matter, occipital white matter).
For automated measurement of MS lesion susceptibility value, the bright lesion signal on T2 weighted (T2w) images can be used to segment lesion. White matter can be automatically segmented on T1 weighted images. A seed point can be placed automatically accordingly the lesion center of mass on T2w. An iterative region growing and thresholding (IRGT) algorithm can be used to segment the lesion volume, which is at the iteration step when the IRGT volume stops growing as the growth has traced the whole lesion and then starts to jump to a large value as the growth leaks to a surrounding structure.
We examined MR images of MS patients from August 2011 to January 2015 with at least two successive MRI sessions that included T2-weighted (T2w), Gd-enhanced T1-weighted (T1w+Gd) and GRE imaging. QSM was constructed in an automated manner from GRE data by deconvolving phase with the dipole kernel that connects tissue susceptibility with the magnetic field estimated from MRI phase. We compared the lesions on two successive MRIs and identified patients with at least one new T2w lesion, i.e., a lesion that was not present in prior brain MRI (in a follow-up MRI that was <1 year from the baseline MRI). All the new lesions were then grouped into enhancing and non-enhancing on T1w+Gd images.
All examinations were performed on a 3.0-T MR scanner (Signa HDxt; GE Healthcare, Milwaukee, Wis.), with an eight-channel head coil. The sequences for each patient were: (a) T2w fast spin echo, (b) pre- and (c) post-gadolinium 3D inversion recovery-prepared Tiw fast spoiled gradient echo, (d) three-dimensional T2*w spoiled multi-echo GRE sequence. Imaging parameters for multi-echo GRE sequence were as follows: repetition time, 57 msec; number of echoes, 11; first echo time, 4.3 msec; echo time spacing, 4.8 msec; flip angle, 20°; bandwidth, 244 kHz; field of view, 24 cm; matrix, 416×320; section thickness, 2 mm. The GRE sequence was performed before Gd injection. The total imaging time was 16 minutes 30 seconds.
QSM was constructed from GRE data using the morphology-enabled dipole inversion (MEDI). The images obtained by the other modalities were registered to QSM by using FSL software (FLIRT Linear Registration Tool, FMRIM Software Library, Oxford, England).
After localizing all new T2w lesions by comparing them with their previous MRIs, three neuroradiologists (J.C., A.G. and G.C, with 18, 9 and 8 years of experience, respectively) used the T1w+Gd images to classify those lesions as enhancing or non-enhancing. They also classified all lesions on QSM as hyperintense and isointense relative to the adjacent white matter. All differences in lesion classification were resolved by the majority.
One neuroradiologist (Y.Z., 4 years of experience) drew the areas of each localized lesion on the T2w images while blinded to the Gd enhancement classification. White matter regions without abnormal signal on T1w and T2w images were identified as normal appearing white matter (NAWM). For a zero reference, a region-of-interest (ROI) was chosen on the NAWM at the contralateral mirror site of an identified lesion with a similar shape and size on T2w images. Then, the ROIs of lesions and NAWM references were overlaid on the QSM images using a semi-automatic software to assess values of lesion susceptibility. Veins or artifacts inside the ROIs were excluded by inspection.
Using relative susceptibilities as a means for distinguishing enhancing from non-enhancing lesions, receiver operating characteristics (ROC) were assessed to determine sensitivity, specificity, and the optimal cutoff susceptibility value (in parts per billion (ppb)). Bootstrapped estimates of the AUC and 95% confidence intervals were produced to evaluate variance. The jackknife cross validation technique was used to evaluate predictive performance of the model. A generalized estimating equation (GEE) was used to predict QSM values from three lesion types: nodular, shell, and non-enhancing. This model assumes a Gaussian distribution and an exchangeable correlation structure to account for the multiple lesions per patient. The GEE analysis was also used to predict QSM values from enhancing and non-enhancing lesions, accounting for repeated measurements per patient. All statistical analyses were performed using SPSS for Windows (version 16.0; SPSS, Chicago). P<0.05 was considered statistically significant. The accuracy for identifying patients with enhancing lesions was also calculated.
From the eligible 482 MS patients, we identified 54 patients with at least one new T2w lesion; there was a total of 133 new T2w lesions. (One patient was excluded because of motion artifacts on GRE images.) The mean age of the 54 remaining patients (12 men and 40 women) was 34.7 years±8.1 (range 20-52 years). The disease duration for these patients ranged from 0 to 18 years (5.71 years±4.51) and the expanded disability status scale scores (EDSS) ranged from 0 to 6. Table 1 shows the demographics of these patients.
On T1w+Gd, 86 (64.7%) of the 133 lesions from 33 patients were identified as enhancing, and 47 (35.3%) non-enhancing from 25 patients, with complete agreement among the three readers. For enhancing lesions, 69 (80.2%) out of 86 were found to be isointense on QSM and 17 (19.8%) slightly hyperintense in contrast to adjacent white matter. According to their enhancement on T1w+Gd, the enhancing lesions were divided into 69 nodular and 17 shell. 13 of the 17 hyperintense enhancing lesions were shell enhancing. All 47 non-enhancing lesions were hyperintense on QSM but 4 (8.5%) of them were only slightly hyperintense. Example images are illustrated in
The lesions' mean susceptibility relative to NAWM was 20.26±7.55 ppb (mean±SD) for non-enhancing lesions and 2.49±6.39 ppb for enhancing lesions (both nodular and shell). In the GEE analysis of lesion susceptibility values between the three lesion types, both nodular enhancing (beta=−19.6, 95% CI=−23.5 to −15.8, p=<0.0001) and shell enhancing lesions (beta=−13.5, 95% CI=−19.0 to −8.0, p=<0.0001) had significantly lower susceptibility values as compared to those of non-enhancing lesions. In the GEE analysis of susceptibility values between enhancing and non-enhancing lesions, enhancing lesions had significantly lower susceptibility values as compared to those of non-enhancing lesions (beta=−17.2, 95% CI=−20.2 to −11.2, p=<0.0001). The exchangeable correlation coefficient was 0.12 for the lesion susceptibility model.
The ROC curve constructed from mean relative susceptibility values of lesions is shown in
QSM0 method is proposed using automated ventricle segmentation and adding L2-regularization of ventricular CSF susceptibility variation to current single orientation MEDI QSM. Our data demonstrate that the QSM0 method reduces artifacts, including CSF susceptibility inhomogeneity and hypointense shadow, which may be induced by susceptibility anisotropy in MEDI. Meanwhile, QSM0 does not affect MEDI tissue contrasts such as MS lesions relative to NAWM. Therefore, QSM0 enables a consistent CSF zero reference.
The zero reference is needed in QSM for absolute susceptibility values, because the measured tissue field is invariant with any uniform shift in the susceptibility distribution. Because CSF is chemically similar to pure water, it is a common choice for the zero reference, using a manually drawn ROI in the atrium of the lateral ventricles. However, the size of the ventricles strongly varies across subjects, and CSF often appears non-uniform in QSM. Consequently, the measured CSF value is sensitive to the ROI choice, bringing an inter-subject uncertainty in the zero reference. The described QSM0 method addresses this challenge in two aspects: First, it eliminates the need for manually choosing CSF ROI, since the ventricles are automatically segmented and used as zero reference regions; Second, it enforces CSF homogeneity during reconstruction, thereby decreasing the sensitivity to the ROI, as is evidenced by the smaller standard deviations within CSF.
The insight for the described CSF-specific regularization comes from studying the origin of the CSF susceptibility inhomogeneity in current QSM, which may be induced by the flow through lateral ventricles and anisotropy of the surrounding white matter tracts. In this work, the white matter anisotropy as the cause of CSF inhomogeneity is investigated in a numerical phantom by simulating tissue fields with and without anisotropy effects. The field from an anisotropic component (χ13 and χ23) generates spatially smooth shadow artifacts in QSM, as shown in
The L2-regularization of CSF susceptibility variation is conceptually similar to the use of a weighted Total Variation as a tissue structure prior that helps reduce streaking artifact. The null space of the dipole kernel allows streaking artifacts that are wave propagations associated with granular noise errors in the measured field data. The Bayesian strategy for reducing the streaking is to employ a prior term of penalizing structures (edges) that are not consistent with tissue structures in anatomic (typically magnitude) images. Similarly, the CSF susceptibility inhomogeneity and shadow artifacts may be regarded as artifacts arising from spatially smooth errors in the scalar dipole field model. A corresponding strategy to reduce these artifacts is to employ a prior term that penalizes CSF susceptibility variation, as described in this QSM0 work.
QSM0 may directly improve the delineation of deep gray nuclei, including the subthalamic nucleus, which is an important surgical target of deep brain stimulation. As QSM has already been used in many clinical applications, it is important to examine the extent to which QSM0 alters the relative contrast outside the CSF. For measuring relative susceptibility differences between MS lesions and NAWM (
Our clinical data suggest that QSM and T2w together allow for accurate identification of enhancing lesions in MS patients without Gd injection within new lesions on serial MRI. This may be a potential clinical application of the reported observation that the magnetic susceptibility of an MS lesion increases rapidly as it changes from Gd enhancing to non-enhancing. One application of QSM in MS is that in serial MRI during regular monitoring of MS patients, QSM may substitute for Gd enhancement in assessing inflammatory activity. Enhancement on T1w+Gd is the current standard method to assess ongoing CNS inflammation for monitoring optimizing inflammation-suppressing treatment. Following the initial inflammatory reaction, the BBB opens and immune cells infiltrate the brain for about 3 weeks; therefore, T1w+Gd may only offer a small window into lesion pathology. During this period, the microglia/macrophages (m/M) take up and degrade myelin fragments, this is reflected in the initial lack of change in the susceptibilities of active lesions on QSM. However, after the BBB seals, immune cells remain active in the brain tissue. For example, m/M remove diamagnetic myelin fragments, and at the same time or afterwards, m/M cells with paramagnetic iron gather both at the periphery and within a lesion to further promote inflammation. For these reasons, both myelin debris removal and iron accumulation likely contribute to the increase in lesion susceptibility observed on QSM. MS lesions are hyperintense for a few years, typically with bright rims on QSM, these bright rims can be interpreted as iron. Therefore, another QSM application in MS is that, including QSM rather than Gd enhancement alone, in an MRI protocol for MS patients may provide more detailed insight into early lesion dynamics in MS.
In the proposed automated segmentation process, we selected the parameter R based on cases in the current study. For cases with different ventricle anatomy such as fetal and pediatric imaging, the parameter might be further tuned for robust segmentation of CSF. The algorithm with effectiveness in reducing “shadow artifacts” in the CSF ventricles can be extended to suppress general shadow artifacts throughout the whole brain region.
Here MLF is an operator to extract the low spatial frequency component in the whole brain region. In one variation of implementing Eq.27, MLF may include a uniform mask of ventricle CSF as in Eq.25, low frequency components of spatial variation at the half of superior-inferior or/and right-left or/and anterior-posterior size of the brain.
In summary, the proposed QSM0 method through automatic segmentation of CSF and enforcement of susceptibility homogeneity within ventricular CSF reduces QSM artifacts from susceptibility anisotropy and provides a consistent zero reference for the QSM value.
In this section, we describe a method to measure liver iron content using QSM to compensating for R2* confounding factors including fibrosis, fat and edema.
The liver is the only organ whose iron content is invariably increased in all types of systemic iron overload, whether from transfusion, increased dietary iron absorption, or both. Measurement of the liver storage iron concentration is the reference method for assessing the magnitude of the body iron burden in the diagnosis and management of both hereditary and acquired forms of iron overload. With iron overload, the amounts of iron in the functional and transport pools changes only slightly; almost all the excess iron is sequestered as paramagnetic ferritin and hemosiderin iron. MRI is highly sensitive to paramagnetic iron and has emerged as the primary non-invasive modality for tissue iron quantification needed in staging iron overload and monitoring iron reduction therapy. The transverse relaxation rate R2 (=1/T2), or the faster and more sensitive R2*(=1/T2*) are readily measured by commercially available MRI pulse sequences to quantify liver iron content (LIC). A major challenge in the R2 and R2* methods for mapping LIC is that they can be confounded by fat, fibrosis and other changes in cellularity that are known to contribute to R2 and R2*. These R2 and R2* confounding effects depend sensitively on the detailed microenvironment of water—confounder interaction and on the confounder spatial distribution within a voxel, which are difficult to model quantitatively. Consequently, R2 and R2* are complex functions of iron and confounders, making it difficult extract the quantity of iron from R2 and R2* measurements. It should be noted that fat suppression does not affect fat contribution to the R2 and R2* of water. In general, R2 and R2* contributions from fibrosis and other changes are difficult to quantify and compensate,
While the magnitude data of the multi-echo gradient echo MRI is used to generate R2*, the phase data can be used to generate quantitative susceptibility maps (QSM) to directly measure magnetic susceptibility without blooming artifacts in R2*. QSM has recently been applied to abdominal organs by additionally accounting for the fat contribution to the measured signal phase. Magnetic susceptibility is an intrinsic tissue property that is directly linearly related to liver iron. The fat contribution to tissue susceptibility can also be linearly compensated according to chemical composition. Other tissue cellularity changes including fibrosis and edema do not affect susceptibility. Accordingly, QSM is developed here to overcome R2* confounding factors in mapping LIC.
QSM algorithms to date have been designed to minimize the streaking artifacts arising from granular noisy sources using a regularization penalizing edges of streaks that are not in known tissue anatomy, which has been very successful. Following this Bayesian strategy, we propose to minimize shadow artifacts using biological prior information: a uniform region of interest on anatomic images such as T1 and T2 weighted images (MU, which may be automatically identified using the characteristic feature of uniform intensity and can be detected using a regional variance calculation) should have uniform iron concentrations. For example, in thalassemia major patients, cirrhosis is typically absent and liver iron overload is evenly spread through the whole liver. Fibrosis with negligible magnetic susceptibility but with substantial T2 contributions (χ
A to zero. We have obtained preliminary data demonstrating that in brain QSM, minimal susceptibility variation in cerebrospinal fluid in the lateral ventricles can be achieved for a consistent and automated zero reference and for additional image quality improvement. Then hepatic QSM with zero reference and shadow-resistance against motion, flow, noise and other errors can be formulated as:
Me the binary edge mask derived from the magnitude image (105). The fat susceptibility χfat can be estimated from the fat map by using χfat=XOA·ρfat/ρOA with ρfat/ρOA is the fat fraction image with ρOA the mean ρfat value of subcutaneous fat regions that are assumed to be pure oleic acid. In the case of large SNR, e−ib−e−i(d*χ)≈−i(b−d*χ), linearizing the first term as fitting field data to the dipole model. Human adipose tissue is composed largely of oleic acid with a volume susceptibility XOA=0.75 ppm. The pure fat may automatically be identified using maximal signal regions on the fat image. The fat contribution to susceptibility can be compensated according to fat fraction and pure fat volume susceptibility to generate quantitative iron map. Absolute LIC for liver iron biomaterial can then be obtained from fat-corrected susceptibility χc using the ferritin/hemosiderin specific mass susceptibility X=1.6·10−6 m3/KgFe:LIC=χc/X.
Once the liver iron concentration is estimated from the fat corrected QSM as described in the above, then liver iron contribution to R2* can be estimated according to calibration curves C(χc) for various imaging parameters and iron concentration. R2* increase from surrounding healthy liver tissue ΔR2* comes from iron and fibrosis. Therefore, liver fibrosis biomaterial can be quantitatively mapped from R2* with an iron correction: ΔR2*|fib=ΔR2*−C(χc).
The L2 biological constraints in Eq.28 may greatly improve QSM reconstruction SNR, as we have found in the piece-wise constant model. Accordingly, this hepatic QSM based on liver morphological priors should allow robust and accurate LIC measurement.
To validate Gd quantification in the presence of fat, four fat-water-gadolinium (Gd) phantoms were constructed comprising of twenty-four balloons with a 28-mm diameter in varying fat/water concentrations with proton-density fat-fraction (PDFF)=0, 14.3, 28.9, and 43.7% and Gd concentration=1.25, 2.5, 5.0, 7.5, and 10.0 mmol/L. Each fat-water-Gd solution included mayonnaise (Olinesa Ltd., Bulgaria) and Gd solutions (Shanghai Xudong Haipu Pharmaceutical Co., Ltd, China). Gd solutions with Gd concentrations of 1.25, 2.5, 5.0, 7.5, and 10.0 mmol/L are corresponding to susceptibilities of 0.41, 0.81, 1.63, 2.45, and 3.26 ppm at room temperature, assuming a molar susceptibility of 326 ppm L/mol at 293 K. The balloons were embedded in a plastic cuboid container filled with water.
Nine controls with no known hepatic disease, nine patients with known or suspected hepatic iron overload, and nine patients with focal liver lesions participated in this IRB approved HIPAA compliant study. All participants provided written informed consent. Subjects with iron overload had either hereditary hemochromatosis, transfusional hemosiderosis, leukemia, or a myelodysplastic syndrome.
The fat-water-Gd phantoms were studied using a clinical 3.0T MR imaging system (MAGNETOM Trio Tim; Siemens Healthcare, Erlangen, Germany) with a combination of head and neck coils. QSM images and R2* mapping of the phantom were obtained from the same 3D unipolar readout multi-echo gradient-echo sequence with the following parameters: TR=13.1 ms, TE1=1.24 ms, ΔTE=2.08 ms, Number of echoes=6, receiver bandwidth (RBW)=1028 Hz/pixel, FA=3°, Matrix=256×208×48, spatial resolution=1.5 mm×1.5 mm×2.0 mm.
The subjects with no known hepatic disease, and those with known or suspected hepatic iron overload were studied using a clinical 1.5T MR imaging system (MAGNETOM Aera; Siemens Healthcare, Erlangen, Germany) with a combination of spine and torso coils. QSM images and R2* mapping of the liver were obtained from the same three dimensional (3D) bipolar readout multi-echo gradient-echo sequence with the following parameters: TR=10.0 ms, TE1=1.44 ms, ΔTE=1.36 ms, Number of echoes=6, RBW=1065 Hz/pixel, FA=6°, Matrix=224×168×44, spatial resolution=1.7 mm×1.7 mm×5.0 mm. Additionally, a generalized auto-calibrating partially parallel acquisition with acceleration factor of 2 in the anterior-posterior direction and elliptical sampling were used to reduce acquisition time. The scan required 14 seconds that was readily completed during a single breath hold.
The subjects with focal lesions were imaged using a clinical 3.0T MR imaging system (MAGNETOM Verio; Siemens Healthcare, Erlangen, Germany). QSM images and R2* mapping of the liver were obtained from the same 3D unipolar readout multi-echo gradient-echo sequence with the following parameters: TR=14.0 ms, TE1=1.49 ms, ΔTE=2.06 ms, Number of echoes=6, RBW=1502 Hz/pixel, FA=10°, Matrix=256×176×26, spatial resolution=1.5 mm×1.5 mm×5.0 mm. The scans required 20 seconds and were completed during a single breath hold.
Phase images were extracted from the complex MRI data. Phase discrepancies between odd and even echoes in the bipolar readout gradients caused by non-ideal gradient behaviors were first measured, modeled as polynomials in space, and corrected in field mapping. Then the corrected complex MRI signal s({right arrow over (r)}, t) was modeled according to tissue chemical composition of water and fat: water signal
with water fraction ρw in a voxel and fat signal
with fat fraction ρf in the voxel, in addition to the phase factor e−i2πb({right arrow over (r)})t corresponding to the magnetic field generated by susceptibility. Therefore,
Simultaneous phase unwrapping and chemical shift determination using graph cuts with conditional jump moves was then performed on the phase images to generate an estimate for fat and water fractions and the magnetic field, which were used as initialization for iterative computation in water-fat signal modeling to fine-tune the field map. The output field was then processed with background field removal using the projection-onto-dipole-field method, and the remaining susceptibility field was processed to generate QSM using the morphology enabled dipole inversion algorithm. The output field was also processed using the preconditioned total field inversion without background field removal.
For an effective R2* mapping of the liver, a Levenberg-Marquardt non-linear fitting was used to fit the magnitude component of the multi-echo data with multipeak fat.
Liver magnetic susceptibility was measured on QSM using region of interest (ROI) analysis. The latissimus dorsi muscle was the reference for zero susceptibility. ROIs of approximately 12 cm2 were placed on the liver region away from major vessels. ROIs of approximately 4 cm2 were placed on the nearby latissimus dorsi muscle. ROIs for focal lesions were placed in the lesions, with an area corresponding to just below lesion size.
The fat contribution to liver susceptibility was estimated and removed from the measured liver susceptibility in the following manner to generate fat corrected QSM. Intensity inhomogeneity due to receiver coils was first estimated using a level set method and was then used to correct the fat image.
Statistical analyses were carried out using SPSS for Windows version 17.0. Linear regression analysis was performed to correlate liver magnetic susceptibility with liver R2*.
Water/fat separation was successful in all 9 controls and 17 of 18 patients. It failed in one patient with extremely high hepatic iron overload. The GRE images from this patient had a detectable liver signal only in the first echo and were therefore excluded from further analysis. Focal lesions in patients consisted of three hemangiomas, one hepatocellular carcinoma (HCC), three suspected cholangiocarcinomas or metastases, one with multiple malignant tumors, and one with multiple liver metastases after interventional therapy. The focal lesions had substantial effects on R2* but minimal effects on fat-corrected QSM. Examples are detailed in
In the patients with focal lesions, there was a correlation between liver R2* measurements and fat percentages (
For ROI measurements that avoided focal lesions and vascular structures, liver susceptibility and R2* measurements were linearly correlated (R2=0.797, P<0.0001,
These data on liver iron mapping in 27 subjects demonstrates that the liver pathology including cysts, fibrosis, hemangioma, hepatocellular carcinoma, cholangiocarcinoma and fatty infiltration confounds R2* measurements but does not affect QSM. Lesions without fat or hemorrhage may be decreased on R2* but blend in with background liver on QSM. The fat contribution to susceptibility was corrected on QSM according to a linear tissue composition model. Our results suggest that QSM can overcome confounding factors in R2* for quantifying liver iron content (LIC).
A major challenge in current MRI measurement of liver iron content is to correct for R2 and R2* confounding factors, including fat, fibrosis and other cellular changes that are often associated with tissue damage from iron overload. While fat fraction can be well quantified in the presence of iron using primarily the unique chemical shifts of fat, fat confounding effects on R2 and R2* depend on the detailed microenvironment of water—fat interaction and the fat spatial distribution within a voxel, which are difficult to model (
Fortunately, MRI signal phase allows rigorous biophysical modeling for the determination of tissue magnetic susceptibility, which is QSM. As a molecular electron cloud property linearly proportional to molecular concentration, the standard linear chemical decomposition can be applied to tissue susceptibility, enabling iron quantification with direct biophysical definition that is not possible with R2 or R2* models of iron quantity. This situation is similar to quantification of deoxygenated heme for mapping cerebral metabolic rate of oxygen consumption (CMRO2), where QSM vastly simplifies the poor traditional empirical nonlinear R2* model. Similarly, QSM enables correction for susceptibility contributions from other cellularity sources to accurately measure the amount of iron. This correction task has been made simple by the fact that fibrosis and most other confounders have virtually zero susceptibility relative to water.
Consistent with the literature, we observed decreased R2* in the presence of liver metastases, hemangioma (
Liver tissue susceptibility in this work was modeled from paramagnetic iron (gadolinium substitute in the phantom) and fat. Through signal phase analysis, the chemical shift of fat allows fat quantification through iterative water fat signal modelling and subsequent QSM. The intravoxel mixture of fat and water with different susceptibilities causes additional field inhomogeneity within a voxel and increases the apparent R2* that depends on imaging parameters including voxel size, field strength and echo time. The strong fat contribution to R2* may explain the correlation between the lesional R2* measurements and the fat percentages (
The preliminary implementation can be improved with the following considerations. First, the number of patients in the preliminary study was relatively small. Patients with location lesions in this study were useful in validating QSM's effectiveness for overcoming R2* confounding factors in mapping the liver iron content. Patients with diffuse liver fibrosis caused by diffuse iron overload are of great clinical concern and should the recruitment focus in future studies. Second, a multi-peak fat model was not considered for QSM. Multi-peak IDEAL (iterative decomposition with echo altering and least-squares-fitting) for water fat decomposition may be adopted in the last fine-tuning step. In one embodiment, the iterative estimation of fat and water fractions and the magnetic field can be initialized using simultaneous phase unwrapping and chemical shift determination based on graph cuts. Another initialization approach that can be used is to include in data acquisition in-phase echoes (the Dixon method for a single-peak model) that allow estimation of R2* and magnetic field for initialization. Third, water/fat separation failed in a patient with very high hepatic iron overload and consequent extremely rapid signal decay. Ultrashort TE techniques based on radial acquisitions or offset echo GRE sequences can be used to acquire sufficient signal for data processing in patients with very rapid T2* signal decay. Fourth, there is no correlation with immunohistochemistry, which would help us interpret QSM and R2* similar to studies on white matter lesions in multiple sclerosis. Fifth, there seemed to be complications resulting in nonlinear susceptibility at high concentrations of Gd and fat in the phantom in
In summary, R2* mapping of liver iron content is confounded by fat, fibrosis, edema and other cellularity pathologies. QSM overcomes confounding factors in R2* for mapping liver iron content.
In this section, we describe a QSM method to measure bone mineralization by detailed biophysical modelling of fat and cortical bone signal.
Magnetic susceptibility is a fundamental tissue property which can be observed in MRI. Densely calcified tissues such as bone have strong diamagnetic susceptibility. Quantitative susceptibility mapping (QSM) can provide quantitative, reproducible images of magnetic susceptibility sources for assessing the health and disease of most tissues but has been limited in bones. Given the importance of bone mineral density for assessing bone fracture risks in postmenopausal women and the elderly population, QSM may become a useful diagnostic tool for non-invasive imaging of bone health without the use of ionizing radiation.
QSM of the bone has been challenging because it requires complete measurements of phase everywhere within the region of interest (ROI) while cortical bone typically has very low signal at conventional echo times in gradient echo (GRE) imaging. Although water is abundant in cortical bone (˜15% by volume), it mostly exists in “bound” form, connected to the crystalline mineral structures or the collagen matrix. As a result, bound bone water has an ultra-short apparent transverse relaxation time (T2*˜300 μs), resulting in no meaningful phase for QSM reconstruction on conventional MRI. Due to these limitations, previous work in muscular-skeletal applications of QSM was either focused on cartilage, or utilized piece-wise estimations of bone susceptibility. An additional problem is the intermingling of fat and water protons in the bone marrow.
In this work, we measured bone MRI signal using an ultra-short echo time (UTE) pulse sequence with TE<<1 ms. We investigated chemical shift and R2* components to properly model bone MRI signal in both UTE and conventional GRE. Finally, we present an MRI method for quantitatively mapping bone susceptibility from UTE-GRE datasets.
Magnetic field estimation is an essential initial step for QSM. For QSM of anatomical regions where magnetic susceptibility is the predominant contributor to the proton phase accrual such as in the brain, field can be estimated from the MRI signal phase. For QSM on regions with high fat content, such as bone with marrow, a complex MRI signal is modelled according to tissue chemical composition of water and fat. The inhomogeneous field generated by susceptibility is estimated through robust separation of water and fat signals.
In a multi-echo GRE sequence, the temporal behavior of the signal originating from the voxel {right arrow over (r)} that contains multiple species can be expressed in general form as:
Here, ρk is the complex amplitude of the signal originating from the kth chemical species within the voxel at time t=0, αnk is the relative amplitude of the nth peak in the chemical spectrum of the kth species, ƒnk and R*2
In many clinical applications, water and fat are the dominant species and Eq.29 becomes
Here αnf is the relative amplitude of the nth spectral peak of fat, ρf and ρw are fat and water complex amplitude at time t=0, R2
Assuming all fat spectral peaks experience similar R2* decay, Eq.30 can be simplified as
Currently, the difference in water and fat R2* within a voxel {right arrow over (r)} are commonly neglected:
Unfortunately, for bone tissue, water is distributed through porous mineral matrix, and fat in marrow tends to aggregate together, and we have R2
Therefore, we propose to avoid the R2* simplification in Eq.32. The water signal is the main component acquired in cortical bone voxels. In our experiments (see below), the largest echo time is much shorter than the T2* of surrounding soft tissues, including marrow. Therefore, the bone water is the only species experiencing significant decay during acquisition. Accordingly, we further propose a bone-specific R2* model with the following reduction of parameters:
Pulse sequence. A radial 3D GRE UTE sequence with TE≥40 μs was implemented on a clinical 3T scanner (GE Excite HD, Milwaukee, Wis.). The sequence used a nonselective hard pulse (rectangular pulse of 100 μs duration) to achieve volumetric excitation and two readouts per TR to accelerate acquisition (
Porcine phantom experiment. A phantom consisting of a porcine hoof (length 16 cm, thickness 6.5 cm) embedded in 1% agarose gel (container height 20 cm, average diameter 11 cm) was constructed and for the purpose of comparison was imaged using an MR scanner (3T, GE Excite HD, Milwaukee, Wis.) and a clinical CT system (LightSpeed Xtra, GE Milwaukee, Wis.). MR imaging parameters included the following: TE=0.04, 0.24, 3.0, 4.0 ms, TR=12 ms, FA=15°, FOV=18 cm, 32000 radial projections per echo, voxel size 0.7×0.7×1.4 mm, BW=±62.5 kHz, total acquisition time=13 min. CT data was acquired using the following parameters: 120 kVp, 200 mA, 0.625 mm slice thickness, 512×512 matrix, achieving nearly isotropic resolution.
Volunteer MR imaging. After obtaining written informed consent, 7 healthy volunteers (6 male, 1 female age range 25-32 years, mean=30) were scanned. All exams were performed under protocols approved by the Institutional Review Board and were HIPAA compliant. Images (11 distal femurs, 1 distal radius and 1 proximal tibia) were obtained on a 3T MRI system using an 8 channel transmit-receive knee coil and 8 channel transmit-receive wrist array. The 3D UTE GRE sequence was used with the following parameters: TE=0.04, 0.24, 3.0, 4.0 ms, TR=12 ms (pairs of the 1st/3rd and 2nd/4th echoes were acquired on successive TRs in this experiment), FA=15°, FOV=18 cm, 32000 radial projections, voxel size=0.7×0.7×1.4 mm3, Total acquisition time for all echoes=13 min.
Image processing and analysis. Each radial MR dataset was reconstructed using regridding, which interpolates the measured signal from radial spoke onto a Cartesian grid. In this work, NUFFT with Kaiser-Bessel kernel and min-max interpolation was implemented. Because of the non-uniform sampling density of our radial trajectory, density compensation was applied to the measured signal prior to re-gridding. Each data set was reconstructed on a coil-by-coil and echo-by echo basis.
For each coil, the phase of the 1st echo was subtracted from the phase of each following echo in order to remove the phase of the coil sensitivity. Then for each echo, the complex images after the phase subtraction were summed across all coils. Finally, the phase of this summed signal was used as the phase for the corresponding echo. Phase of the 1st echo became zero after this process. In mathematical form this process can be written:
Î
m
={circumflex over (M)}
mexp(−i{circumflex over (Φ)}m)
where Îm is complex data corresponding to the mth echo, {circumflex over (M)}m=√{square root over (|Σn=1N
Water, fat and inhomogeneity components were then obtained using IDEAL-based techniques. For this iterative method, it is important to use an initial guess for the inhomogeneity field that is reasonably close to the true solution in order to avoid convergence to a local minimum. In the current work, a preliminary field map estimation was carried out using a graph cut based phase unwrapping and chemical shift estimation assuming a single chemical shift peak with value ƒ=−3.5 ppm·γB0 Hz. An initial estimation of the distribution of R2* values was produced by mono-exponential fitting. The complex data was then fitted to Eqs [33] and [32] assuming different models for the fat chemical spectrum—single peak model (Δƒ=−3.46 ppm) and multi-peak model ({right arrow over (Δƒ)}={−3.82, −3.46, −2.74, −1.86, −0.5, 0.5} ppm). The resulting inhomogeneity field b was then passed for susceptibility mapping to the morphology enabled dipole inversion (MEDI) pipeline in which Projection onto Dipole Fields (PDF) was utilized for background field removal.
The CT images of the porcine hoof were resampled and registered to the reconstructed susceptibility map using the FLIRT algorithm in the FSL toolbox. ROI analysis was performed on co-registered volumes to correlate the CT signal (in Hounsfield units) with the MR signal (calculated susceptibility values). For this, 58 in-slice regions of interest were manually drawn in trabecular and cortical areas of metacarpal and phalanx bones, and flexor tendon. QSM (referenced to muscle tissue) and CT volumetric averages for the ROIs were recorded and used for correlation analysis.
Porcine phantom results. A susceptibility map of the specimen was successfully reconstructed. Comparison of CT and QSM images is shown in
Volunteer results. QSM reconstruction using the proposed method with water-specific R2* and multi-peak fat spectrum modelling was successful in all volunteers (see example in
The results of the volunteer scans (Table 1) reveal good intra- and inter-subject agreement of susceptibility values for cortical bone (inter-subject average χbone=−1.4±0.2 ppm).
Our data demonstrate the feasibility of mapping bone magnetic susceptibility including cortical bone, trabecular bone, and marrow. Bone QSM is improved by utilizing a radial UTE gradient echo (GRE) sequence to acquire the rapidly decaying signal of bone water, and a novel bone-specific-R2* in signal modeling to achieve accurate water/fat separation. Results in an ex vivo porcine hoof show excellent correlation with CT. By comparing QSM reconstructions using the current standard signal model with a common R2* for both fat and water components, the bone-specific-R2* approach reduces fat fraction estimation errors within cortical bone and generates realistic susceptibility values for cortical bone in 7 healthy volunteers with excellent depiction of all skeletal tissue components.
Field estimation from UTE and conventional echo GRE data using the conventional IDEAL is non-trivial, prone to noise propagation, and R2* bias. The robust extraction of biophysical parameters from GRE data require proper signal modeling with a minimal number of parameters. Standard GRE signal models assign identical R2* to both water and fat species fail to correctly estimate the field, consequently generating bone QSM with grossly erroneous values (such as paramagnetic bone), because rapid signal decay of bone is misinterpreted as dephasing due to the high fat fraction. We realize that over the echo times acquired in this study, there is negligible R2* decay of the fat signal but rapid R2* decay of water dispersed within the calcified matrix. This leads us to the bone-specific R2* model that accurately reflects the underlying tissue physics, attributing signal decay exclusively to the water component. Both the porcine hoof and human volunteer data show that this bone-specific R2* model when combined with iterative least squares fitting provides reliable estimation of chemical species distribution as well as a biologically plausible field map. Additionally, the application of the multi-peak fat spectrum model further improves quality of the water/fat separation and, consequently, the susceptibility map.
Previous studies show bone magnetic susceptibility can be estimated by a highly regularized piece-wise inversion of the local magnetic field, which may be problematic when there is a local variation in tissue susceptibility. A masked iterative calculation of susceptibility can also be used to estimate distribution of strong magnetic sources in tissue, but the iteration tends to enforce both uniformity and underestimation of high susceptibility structures. Above all, these attempts do not properly account for the underlying fat water biophysics in the MRI signal generation, nor do they acquire signal within bone using UTE, as illustrated in this work. For trabecular bones with densities much lower than those of cortical bones, QSM may be straightforwardly applied using conventional echo times.
Bone mineral density (BMD) assessment is central to the diagnosis of osteoporosis. Currently, BMD is measured using dual energy X-ray absorptiometry (DXA), which may be confounded by factors including degenerative changes, small bone size, and overlaying anatomic structures. Quantitative CT (QCT) has been proposed to overcome these limitations of DXA, but it requires a much higher X-ray dose. A number of non-invasive non-ionizing radiation based methods for characterizing bone tissue are also available. Quantitative ultrasound bone characterization is an important modality for clinical assessment of osteoporotic patients. Several studies have reported use of velocities of axial propagation of ultrasound waves along the shaft of long bones as a biomarker for monitoring of bone fracture healing and bone assessment in osteoporotic patients. MRI-based techniques (such as ZTE and UTE) rely on imaging water bound to the organic matrix of the bone. These acquired signals can provide information on density and porosity of the bone material; this requires separation of bound and free water signal components that have significantly different relaxation constants. UTE imaging with bicomponent analysis has been successfully applied before for ex- and in-vivo cortical bone imaging.
Conventional MRI may be used to assess bone quality but not mineral density. The MRI-based bone QSM presented here may provide BMD assessment without the X-ray radiation. Therefore, our results warrant future research in comparing bone QSM with DXA on patients—including those suffering from osteoporosis—to establish an accurate non-X-ray BMD assessment for predicting fracture risk and guiding therapy. Since comprehensive comparison of QSM as a biomarker with other non-ionizing radiation based measures falls outside of the scope of this paper, additional investigation is required to further assess clinical utility of bone QSM.
The preliminary implementation of including UTE in QSM can be improved with the following consideration. The scan time may be reduced substantially using data acquisition acceleration strategies in Bayesian MRI. This preliminary implementation also assumed that water is the only species experiencing significant decay by the time the last echo is acquired. Although correct in many practical situations when bone matrix bound water is present, this assumption may break down under conditions of strong field inhomogeneity (e.g., imaging near abdominal cavity). These field variations can lead to strong spin dephasing and, as a result, rapid R2* decay in soft tissue regions. In this situation, high order shimming might be required. Additional limiting factor not taken in the account in the present work is the image distortion due to off-resonance frequency of fat. Due to the limited bandwidth, in the radial acquisition used in this study off-resonance manifests itself in blurring of the water/fat tissue interfaces in both magnitude and phase images. This blurring is expected to affect final QSM results, and application of off-resonance correction methods appropriate for non-Cartesian acquisition may be required. Finally, as has previously been reported, QSM is prone to underestimation of strong susceptibility sources, and further adjustments of the imaging protocol (matrix size, resolution) may be necessary to minimize error.
Comparison of the results obtained in the current study with those previously reported in the literature demonstrates possible underestimation of bone magnetic susceptibility obtained with the proposed technique. For example, in authors report bone susceptibility χbone≈−2.4 ppm obtained in in vitro experiment. Further, in vivo estimation performed in yielded fairly wide range of values χbone∈{−1.8, −2.3} ppm depending on the used reconstruction method. Most notably, this underestimation can be observed around the joints (
The bone susceptibility, being proportional to CT attenuation coefficients, can be directly used to correct for attenuation in positron emission tomography reconstruction. This would be especially useful for the combined MRI and PET scanner (MR/PET).
In conclusion, quantitative magnetic susceptibility maps across the entire bone cross-section are feasible using a combination of UTE, a conventional TE gradient echo acquisition, and a bone-specific R2* signal model.
Implementations of the abovementioned techniques can be used in a variety of applications. As examples, the following section describes the application of one or more of these techniques to numerical simulation, phantom experiments, and human brain imaging.
Aerobic respiration is the main energy supply for the brain: it is ˜2% of the adult human body weight, yet it consumes 20% of total oxygen supply. Consequently, oxygen deficiency can easily damage brain tissue, as in hypoxia in Alzheimer's disease (AD) and multiple sclerosis, and ischemia in stroke. qmCMRO2 is specifically valuable for evaluating these brain disorders. For example, in advanced stroke therapy, it is key to identify ischemic penumbra that can be salvaged. Cerebral blood flow (CBF) reduction is problematic for defining the ischemic core. Diffusion perfusion mismatch fails to demonstrate improved clinical outcomes. Oxygen extraction fraction (OEF) elevation in the peri-infarct region is only a weak indicator of tissue outcome. These inadequacies are overcome by CMRO2(=CBF*OEF*arterial oxygen), which consistently predicts neuron death. Another example is the altered aerobic glycolysis in AD. The measurement of aerobic glycolysis depends in part on the quantitative measurement of oxygen metabolism. The growing concerns of AD in the aging population has generated increased interest in qmCMRO2.
In addition to examining disease, qmCMRO2 has also been an instrumental tool for studying basic brain function. For example, the discovery of the default mode network is arguably the most important advance in our knowledge of brain function in the last two decades, and it has relied on the use of qmCMRO2.
The 15O tracer based PET imaging technique has been pioneered by researchers at Washington University (WU) and remains the gold standard for quantitative imaging of brain circulation and metabolism. PET using 15O-oxygen, 15O-carbon monoxide, 15O-water and 18F-fluorodeoxyglucose (FDG) provides the only truly quantitative approach to measure human brain circulation and metabolism in health and disease. Although blood oxygenation level dependent (BOLD) fMRI is widely used for qualitative mapping of brain function, the quantitative interpretation of such mapping is difficult. Unfortunately, the short 2-min half-life of 15O necessitates onsite production and limits its dissemination in both clinical and research settings.
To address the limitations of 15O PET, MRI has been investigated for qmCMRO2. CBF can be mapped using spatially selective RFs to generate arterial spin labeling (ASL) or using a contrast agent to label spins, but qmCMRO2 requires additional estimation of deoxyheme concentration ([dH]) from MRI. So far, efforts have been focused on magnitude signal that unfortunately has a complex dependence on [dH]. One insightful idea is to estimate [dH] by mapping blood T2 in the venules using velocity selective labeling, i.e., the quantitative imaging of extraction of oxygen and tissue consumption (QUIXOTIC) method and its variants. However, this approach suffers from problematic flow capture, poor sensitivity, and contamination from arterial blood and cerebrospinal fluid. The [dH] dependence of MRI magnitude signal has been rigorously modeled for an approximate vasculature, i.e., the quantitative BOLD (qBOLD) method, but there is error in approximation and the inverse condition for estimating [dH] is poor. The most widely investigated approach is calibrated fMRI (cfMRI), which empirically models the gradient echo (GRE) magnitude R2* decay. The [dH] estimation from R2* requires data acquisition during two vascular challenges, typically hypoxia and hypercapnia, which make cfMRI cumbersome to perform in research and difficult to use in clinics. In summary, these MRI based qmCRMO2 techniques have poor sensitivities, use uncertain biological assumptions, or may require complicated vascular challenges. Furthermore, MRI based qmCMRO2 has yet to be validated using PET.
QSM Eliminates Empirical Parameters in Current Calibrated fMRI (cfMRI) Based on R2*.
cfMRI is currently the most widely used MRI approach for qmCMRO2. In cfMRI, the venous deoxyheme concentration [dH] is estimated by modeling the GRE magnitude signal as an R2* exponential decay with R2*([dH]) as:
R
2*([dH])=a·v·[dH]β+r0 [34]
Here, v is the cerebral venous blood volume; r0 is the signal decay rate when [dH]=0. a and β in Eq.34 are empirical parameters without clear physical meaning, as they depend on field strength and vessel geometries, and require additional data acquisitions using vascular challenges such as hypoxia and hypercapnia. Consequently, the cfMRI approach for qmCMRO2 is very difficult to perform in clinical practice.
We can solve this problem by taking advantage of the typically discarded GRE phase signal, which can be used for quantitative susceptibility mapping (QSM). The magnetic susceptibility χ (QSM) model for estimating [dH] eliminates empirical parameters a and β in Eq.34:
χ=v·XdH·[dH]+η0 [35]
Here XdH=151.054 ppb ml/μmol is the dH molar magnetic susceptibility; η0 is the susceptibility when [dH]=0. Therefore, QSM can significantly simplify the cfMRI experiment for qmCMRO2.
Obtaining qmCMRO2 from a QSM analysis of GRE phase data (Eq.35) requires regularization, which can be done using qBOLD with vasculature prior to describe the GRE magnitude data and minimal local variance (MLV) in CMRO2. The uncertainties in the vasculature used in qBOLD and in the MLV of CMRO2 are optimally managed through Bayesian inference. Therefore, qmCMRO2 can be optimally obtained from GRE magnitude and phase data and from ASL data that can be easily acquired on almost every MRI scanner.
Quantitative Mapping of CMRO2 (qmCMRO2) Using QSM cfMRI+MLV+Leenders.
For Bayesian estimation for CMRO2, we approximate both data noise and prior uncertainty as Gaussian. qmCMRO2 requires mapping cerebral blood flow (CBF≡f) which can be obtained using ASL, and the difference between venous dH concentration [dH] and arterial dH concentration [dH]a which is oxygen extraction fraction (OEF≡σ=([dH]−[dH]a)/H) scaled by H=7.53 μmol/ml, the heme molar concentration in tissue blood (hematocrit of Hct=0.3567 (straight sinus)). [dH]a=0.02 H=0.15 μmol/ml for 98% arterial oxygenation. Further, with K=nHbMHbCHbHctXdH=0.0084 L/mol and χ0=Kv[dH]a/H+η0, the CMRO2 equation and Eq.35 become
CMRO2≡ε=Hfσ, χ=Kvσ+χ0. [36]
To optimally estimate CMRO2 ε from ASL data f and GRE phase derived χ, we consider three biological priors: 1) Leenders prior to relate cerebral venous blood volume v to flow data ƒ:v=L(ƒ). 2) In a block much smaller than the Brodmann area or a region of similar cytoarchitecture, gray and white matter may have uniform yet different CMRO2 values. This minimal local variance (MLV) of CMRO2 is consistent with the 15O PET CMRO2 literature and with the fMRI activation maps. 3) The sum of CMRO2 over gray/white matter tissues is equal to the global CMRO2 measurable in large veins. These three priors can be incorporated into the energy minimization problem to estimate ε (and v) as follows:
argminε,v∥w((χ−χ0)ƒ−ΓL(ƒ)ε)∥22+λ1∥v−L(ƒ)∥22+λ2(∈(
)−ε(i))2+λ3∥
Here, wp is the inverse of noise in QSM multiplied by f, Γ=K/H, ∈() is the tissue specific average of ε in a block
, and
)/Hf(i) (block constant ε). Then Eq.37 is simplified to argminε
∥wp((χ−χ0)f−ΓL(ƒ)∈(
))∥22+λ∥
)∥22. By shifting blocks to reduce blockiness, we obtained very encouraging challenge-free MLV results that were comparable with those seen following a caffeine challenge (
The Leenders prior may work well for studying neurophysiology in healthy subjects and may be extended to challenged conditions using the Grubb relation. Disease may decouple blood flow and blood volume, hence, alternative is needed for estimating venous CBV v. Dynamic contrast enhanced MRI only provides total CBV. The vascular space occupancy (VASO) approach can provide v but may be problematic in patients. We will explore the qBOLD approach to over the limitation of the Leenders prior.
Optimal qmCMRO2 Using QSM cfMRI+qBOLD+MLV.
The qBOLD approach uses a prior vasculature of randomly oriented tubes that can approximate veins in a voxel for most healthy and diseased tissues. qBOLD has been applied to generate CMRO2 using multi asymmetric spin echo or GRE. qBOLD only utilizes signal magnitude to estimate multiple parameters: σ and v which are relevant to CMRO2, and S0, R2, χ0 (defined in Eq.38 below) which are not relevant to CMRO2. qBOLD suffers from poor SNR and errors in vasculature assumptions. We propose to utilize the phase data through QSM and the magnitude data through qBOLD. The uncertainties in the qBOLD assumptions and in the MLV and other regularizations can be optimally handled by Bayesian statistics by generalizing Eq.37 as:
argminσ,v,S(σ,v,S0,R2,χ0) [38]
Here, Q(v,σ,S0,R2,t)=S0e−R(v,Gσt)+ve−R
(v,Gσt).
(σ,v,S0,R2,χ0) is the regularization, which includes the MLV in Eq.37, smoothness on S0 and anatomic structural knowledge of R2, χ0. The MLV idea can be extended to the signal histogram space, for example, to impose minimal variance of or similar σ and R2 over all voxels with similar signal behaviors such as a similar temporal decay trend. Eq.38 will be solved using nonlinear Gauss-Newton optimization as in QSM. We will also investigate the use of the primal dual algorithm to speed up the computation.
The study was approved by the local Institutional Review Board. Healthy volunteers were recruited (n=11, 1 female, mean age 34±12 years) for brain MRI on a 3T scanner (HDxt, GE Healthcare) using an 8-channel brain receive coil. All subjects were instructed to avoid caffeine or alcohol intake 24 hours prior to MRI.
MRI was performed before and 25 min after the oral administration of 200 mg caffeine. The protocol consisted of a 3D fast spin echo (FSE) arterial spin labeling (ASL) sequence, a 3D multi-echo spoiled gradient echo (SPGR) sequence, and an inversion prepared T1w SPGR sequence (BRAVO).
The 3D FSE ASL sequence parameters were: 20 cm FOV, 3 mm isotropic resolution, 1500 ms labeling period, 1525 ms post-label delay, 62.5 kHz readout bandwidth, spiral sampling of 8 interleaves with 512 readout points per leaf, 35 axial slices, 10.5 ms TE, 4796 ms TR, 3 signal averages, and 6 min scan time.
The 3D SPGR sequence parameters included: 0.52 mm in-plane resolution, 2 mm slice thickness, identical volume coverage as the 3D FSE ASL sequence, 7 equally spaced echoes, TE1=4.3 ms, echo spacing 7.9 ms, TR 56.6 ms, 62.5 kHz readout bandwidth, 150 flip angle, and 7 min scan time. The pulse sequence was flow-compensated in all three directions.
The inversion prepared SPGR T1w sequence parameters were: 1.2 mm isotropic resolution, identical volume coverage as the 3D FSE ASL sequence, TE 2.92 ms, TR 7.68 ms, prep time 450 ms, ±19.5 kHz readout bandwidth, 150 flip angle, and 2 min scan time.
CBF maps (ml/100 g/min) were generated from the ASL data using the FuncTool software package (GE Healthcare, Waukesha, Wis., USA). For QSM reconstruction, an adaptive quadratic-fit of the GRE phase followed by a projection onto dipole fields was used to obtain the local field map. Susceptibility was computed from the local field map using the Morphology Enabled Dipole Inversion (MEDI) algorithm. Susceptibility values were referenced to the susceptibility of water, as estimated in the cerebrospinal fluid (CSF) masks (see below). All images were co-registered and interpolated to the resolution of the QSM maps using the FSL FLIRT algorithm.
According to mass conservation of oxygen molecules, the cerebral metabolic rate of oxygen CMRO2 (μmol/100 g/min) can be expressed as
CMRO2=CBF·([dH]v−[dH]a) [39]
where CBF is the cerebral blood flow (ml/100 g/min) and [dH]v and [dH]a are the deoxyheme molar concentration (μmol/ml) in venous and arterial blood, respectively.
Oxygen Extraction Fraction (OEF) is the percentage of oxygen extracted from arterial blood and can be expressed as
where [H]=7.53 μmol/ml is the heme molar concentration in tissue blood assuming a tissue blood hematocrit of Hcttissue=0.357, which is estimated from an assumed hematocrit of Hctss=0.47 in the straight sinus and a Hct ratio of 0.759 between large vessel and brain tissue. [dH]a=0.15 μmol/ml assuming 98% arterial oxygenation.
Calculation of heme concentration ([H]) in tissue blood
Where 4 account of 4 heme per Hb, Hcttissue=0.357, ρRBC,Hb=0.34 g/ml is mass concentration of Hb in a red blood cell, MHb=64450×10−6 g/μmol is the molar mass of dHb.
χo=CBV·Xba+CBVa·XdH,mol·[dH]a
Here CBVa=CBV−CBVv=0.23CBV is the arterial blood volume fraction. [dH]a=0.15 μmol/ml is [dH] in arterial blood assuming 98% oxygenation. Xba is susceptibility of 100% oxygenated blood.
X
ba
=HCt
tissueδHb·XoHb+(1−HcttissueδHb)Xp=−108.3 ppb
Where δHb is volume fraction of Hb in RBC.
Here ρRBC,Hb=0.34 g/ml is Hb mass concentration in a red blood cell. ρHb=1.335 g/ml is Hb density in pure aggregate. XoHb is volume susceptibility of oxyhemoglobin (OHb). XoHb=−4πρHb(0.587×10−6 ml/g)−(−9035 ppb)=−813 ppb. 0.587×10−6 ml/g is mass susceptibility of globin relative to vacuum in CGS unit. 4π is for unit conversion from cgs to si units. −9035 ppb is water volume susceptibility relative to vacuum. Xp=−37.7 ppb are the oxyhemoglobin and blood plasma volume susceptibility. [dH]v in Eqs.39 and 40, and subsequently CMRO2 can be mapped using QSM.
Here, χ is the magnetic susceptibility mapped using QSM. χ0 reflects the susceptibility contributions from pure oxygenated blood and arterial deoxyhemoglobin. CBVv is the venous blood volume fraction, estimated as 77% of total CBV. XdH,mol=151.054 ppb ml/μmol is the molar magnetic susceptibility of dH.
In Eq. 41, CMRO2 and χnb are the two unknowns. In previous studies, susceptibility and blood flow of baseline and challenged brain states were mapped. Hence, CMRO2, χnb, and subsequently OEF can be solved. However, CBF challenges such as caffeine administration are time-consuming, presenting difficulties when performing in patients.
To solve this limitation, the QSM and CBF maps were segmented into small blocks and each block was further segmented into two tissue types (gray and white matter) using T1 weighted images. Within each block and each tissue type, χnb and CMRO2 are assumed to be constant. This assumes that, within a small local region, the same tissue type has similar biological and chemical makeup. Hence, CMRO2 and χnb of each block and tissue type can be organized into a vector of unknowns x in a linear system format:
Here ψi=CBVv,i·XdH,mol/CBFi and ϕi=χi−CBVi·(Xba+XdH,mol·[dH]a). The index i=1, . . . , N denotes the ith voxel of a given tissue type with a block with N the total number of such voxels in that block. Xba=−108.3 ppb is the susceptibility of fully oxygenated blood.
In the following, this algorithm will be referred to as minimum local variance (MLV). To reduce error propagation and improve convergence speed, a right preconditioning technique is applied. The new system is defined as APxP=b, where AP=AP and xP=P−1x. P is the preconditioner defined as P=diag(CMRO2,ub, avg(χ)) which is designed to allow the elements of xP to have similar order of magnitude. CMRO2,ub is the upper bound of CMRO2 calculated within each block using Eqs.39&40 using OEF=1. avg(χ) denotes the average of susceptibility of voxels within a block that are of the same tissue type.
The solution for a given tissue type (grey or white matter) was obtained by minimizing the following cost function over the whole brain using a limited-memory Broyden-Fletcher-Goldfarb-Shanno-Bound (L-BFGS-B) algorithm with physiological bounds on CMRO2:
A is diag(AP), a block-wise diagonal matrix with the matrix AP for the given tissue type and for each block along the diagonal, and x, b, and CMRO2 are the concatenated column vectors of xP, b, and CMRO2 for that same tissue type and for each block, respectively. In the second term,
Whole brain grey matter (GM), white matter (WM), and CSF masks were created using the FSL FAST algorithm on the inversion prepared SPGR T1w images. These masks were further divided into 6×6×6 mm3 isotropic blocks. This block based computation was performed for the GM and WM separately, selecting within each block only those pixels that belong to the corresponding tissue type. This allowed a proper convergence when solving Eq.38. To minimize potential bias due to specific block grids, blocks were shifted in (3D) diagonal directions 4 times with a step size of 1.5 mm. The data were reprocessed and the results were averaged to obtain the final images.
CMRO2, OEF, and χnb were also solved as follows. Briefly, assuming χnb remained the same before and after challenge, susceptibility and CBF measurements of each voxel in baseline and challenged states were paired to solve OEF and χnb for every voxels using LBFGS-B solver with identical OEF/CMRO2 bounds and global CMRO2 constraint as the proposed method. The percentage change of CMRO2 between the two states were estimated in the straight sinus and incorporated into the minimization.
For this feasibility study, CBV was estimated from CBF based on a linear regression derived from CBF and CBV measurements in gray and white matter using 15O Positron Emission Tomography (PET) in healthy subjects (n=34): CBV=(0.0723 CBF+1.144)/100. Post caffeine CBV was estimated from baseline CBV using the Grubb's exponent of 0.38 because the empirical CBF and CBV relationship derived from PET data was from healthy subjects without vaso-challenges and may not be suitable for estimating post caffeine CBV. CBV changes after challenge were assumed to be on arterial side as suggested in the literature.
Cortical grey matter (CGM) was segmented into the bilateral anterior (ACA), middle (MCA) and posterior cerebral artery (PCA) vascular territories (VT) by the same neuroradiologist (A.G.). For each ROI, the mean and standard deviation of CMRO2 and OEF values were measured. Bland-Altman plots and paired sample t-tests were performed to compare CMRO2 and OEF values in the VTs using previous and proposed methods. P values lower than 0.05 were considered to indicate statistical significance. All data were expressed as mean±standard deviation.
Numerical simulation was performed to investigate CMRO2 estimation error caused by 1) CBF error, 2) partial volume effect, and 3) QSM error. One block with matrix size of 8×8×5, the same block size in this study, was constructed. The upper and lower half was set to be white (WM) and gray matter (GM), respectively. Inputs were CMRO2, χnb, and CBF (ground truth): CMRO2 of GM: 181, CMRO2 of WM: 138 (mol/100 g/min), χnb of GM: 0.02, χnb of WM: 0.45 WM (ppb). CBF values are given voxel-wise with normally randomized values, 49.5±4.3 for WM and 64.8±3.6 for GM (ml/100 g/min). These values were the average values in a subject and, for the CBF standard deviation, the whole GM and WM standard deviation was calculated and divided by five to obtain a reasonable estimation for the same tissue type in small brain region. The QSM was subsequently estimated voxel-wise with the given CMRO2, χnb, and CBF values. 1) To investigate the CBF error propagation in CMRO2 map, we considered two cases; 1-1) low SNR CBF and 1-2) biased CBF. For low SNR case, noise was added to CBF (SNR: 4.5), and MLV was performed subsequently. This was repeated 100 times to obtain average and standard deviation. CBF SNR was estimated the ratio of mean over standard deviation in the caudate nucleus, 4.5±1.7 (n=5). 2) To investigate the partial volume effect within individual voxels, the QSM and CBF value at the GM boundary layer was averaged with an immediate neighboring WM QSM and CBF value, respectively. This means that GM boundary voxel actually consisted of half GM and the other half WM, but was treated as GM. 3) To estimate QSM error propagation into CMRO2 estimation, noise was added to QSM (SNR 3.1, estimated in the same way as CBF SNR). This was repeated 100 to obtain average and standard deviation. For all cases, global CMRO2 constraint weighting (λ) was set to be 500 as used in this study.
Numerical simulation showed that CMRO2 estimation error was affected by low SNR CBF and by biased CBF. In the case of low SNR CBF (SNR 4.5), CMRO2 error was 6.0% for GM and 4.4% for GM. In the case of biased CBF (10% underestimation), relative CMRO2 error for both GM and WM was 10%. In terms of the partial volume effect at the GM boundary, the relative CMRO2 difference was 3.0% for GM and 2.2×10−5% for WM. In case of the QSM noise propagation into CMRO2 (SNR 3.1), the resultant CMRO2 error was 0.1% for WM 1.5% for GM.
Our preliminary results demonstrate the feasibility of a noninvasive CMRO2 mapping without a vascular challenge based on QSM and CBF using MLV. The CMRO2 measurements agrees well between the challenge and MLV methods in healthy subjects. The elimination of the vascular challenge greatly reduces the complexity of the scanning protocol, the overall exam time and the potential risks associated with the administration of agents that modify blood flow.
Previously reported R2*-based CMRO2 mapping methods are limited by R2* modeling for [dH]. Calibration of the model parameter is required, which may vary with scanning parameters such as field strength. Consequently, measurements at three brain states (normoxia, hypercapnia, and hyperoxia) are needed for quantitative CMRO2 mapping. However, it is very cumbersome to administer stimulants to the brain to generate three distinct states. Magnetic susceptibility has a completely defined linear relation with [dH], eliminating the R2* model parameter calibration and reducing the number of required brain states for CMRO2 mapping to two. The two brain states are needed to separate the susceptibility contribution from blood and non-blood tissue and can be achieved using vascular challenges, such as by administering caffeine or using hyperventilation. However, even administering a single challenge can be cumbersome in practice. Using MLV quantitative CMRO2 map can be obtained without such a challenge, making CMRO2 more feasible to perform in practice.
The MLV implemented in this study used an isotropic block size of 6×6×6 mm3 such that multiple CBF voxels, acquired at a 3×3×3 mm3 isotropic resolution, were available within each block. In 1909, Brodmann divided the human cerebral cortex into 44 areas according to its cytoarchitecture, which were later shown to be related to various neural functions. In this study, the cerebral cortex was divided into more than 4000 regions. Hence, the assumptions of constant CMRO2 and χnb within the same block and tissue type may be considered as reasonable. This MLV implementation then resulted in multiple data points to uniquely estimate the two unknowns CMRO2 and χnb within each block. There is substantial amount of noise propagation from the CBF and QSM input to the CMRO2 output, which can be further reduced using a global CMRO2 constraint as the regularization in Eq.42.
There are several limitations in this study. In MLV, χnb and CMRO2 are assumed to be constant within each block and each tissue type. Therefore, MLV derived CMRO2 values are block-wise constant and depend on accurate tissue segmentations. MLV may be further improved by considering physiologically relevant brain parcellations. The rationale for choosing 6×6×6 mm3 isotropic blocks was that, for this data set, each block contained exactly 2×2×2 CBF and 12×12×3 QSM, both deemed large enough to avoid ill-posedness but not too large to avoid violating the homogeneous tissue assumptions. The block size was chosen to balance between these competing demands. To investigate the effect of block size, we applied different block-sizes (2×2×2 mm3 and 9×9×9 mm3) for one subject. CMRO2 values for cortical GM were 119.0±107.5, 123.5±83.2, and 124.8±64.2 μmol/100 g/min for 2 mm, 6 mm, and 9 mm block size respectively. As the block size increases, the standard deviation of CMRO2 decreases. Averaging over several sub-block shifts, as done in this work, suppresses the resulting piecewise constant appearance. Nevertheless, sharp transitions remain between different tissue zones (
In summary, our study demonstrated the feasibility of a noninvasive CMRO2 mapping without vascular challenge based on QSM and CBF using MLV. In vivo imaging in heathy subjects demonstrates that the CMRO2 and OEF maps obtained with the proposed MLV method agree well with those obtained using a blood flow challenge based method.
In this section, we present QSM of the heart magnetic susceptibility source using navigator gating data acquisition.
Mixed-venous oxygen saturation (SvO2) measures the percentage of oxyhemoglobin that returned to the right ventricular chamber after circulated through the vascular system. It reflects the balance between global oxygen delivery and global oxygen consumption, making it an important measure of cardiopulmonary function that is widely used to manage critically ill patients. Currently, measurement of SvO2 requires an invasive pulmonary artery catheterization procedure. Thus, a method for using non-invasive MRI to measure SvO2 may be of interest for improving patient care.
Quantitative Susceptibility Mapping (QSM) is an MRI contrast mechanism capable of quantifying the difference between diamagnetic and paramagnetic materials. Due to the difference between the number of unpaired electrons in the iron, oxyhemoglobin is diamagnetic and deoxyhemoglobin are paramagnetic. This allows QSM to be used to quantify the oxygen saturation level of blood in the body. However, in vivo QSM has yet to be applied to the heart due to three major challenges: cardiac and respiratory motions, fat chemical shift in the chest and heart, and background field. Here, we will describe our methods for overcoming these challenges to generate in vivo cardiac susceptibility maps, and use it to quantify the relative difference between the oxygen saturation level in the LV blood and RV blood.
Motion has always been a problem for cardiac MRI, and even more so for cardiac QSM, due to the need for long TR. To suppress cardiac motion, we used an ECG triggered sequence, and limited the acquisition window to be less than 200 ms within diastole. To suppress respiratory motion, we asked the volunteers to preform breath-holds during the scan. The volunteers were asked to hold their breaths at maximum exhalation to ensure that each slice was acquired at the same respiratory position.
We also developed a free-breathing navigator method. An ECG-triggered navigator gated multi-echo 3D FGRE sequence was developed for cardiac QSM on a 1.5T scanner (GE Healthcare). A pencil beam navigator echo was used to detect the diaphragm position, and an efficient 2-bin phase-ordered automatic window selection (PAWS) gating algorithm2 was used to control the data acquisition based on the diaphragm position in real time. Two navigator echoes were played each heartbeat, one right before acquisition and one right after acquisition. If the difference between the two diaphragm positions detected by these two navigator echoes was greater than 2 mm, then the data acquired in this heartbeat was rejected. We compared the 2DBH and 3DNAV approaches for cardiac QSM in healthy volunteers.
The scan parameters for the 2DBH approach were as the following: 1st TE/ΔTE/#TE/TR/BW=1.5 ms/2.2 ms/5/15 ms/±62.5 kHz, acquisition resolution=192×192 with 75% phase FOV, 20 slices, 5 mm slice thickness, and 10 views per heartbeat. The scan parameters for 3DNAV were the same except that 24 slices were acquired with 0.75NEX and 0.75 partial slice encoding.
Five healthy volunteers, each consented with IRB approved protocol, were scanned with the 2DBH sequence, in which they were asked to perform breath holds at maximum exhalation, and the 3DNAV sequence, in which they were freely breathing.
The total field was obtained from the complex data from a chemical shift update method3, which was initialized by the output of a graph-cut based phase unwrapping and fat-water separation method4. Then QSM was reconstructed using the Total Field Inversion Method.
A successful dipole inversion requires a continuous 3-dimensional field as the input, so consistent breath holds for every slice is necessary to ensure such conditions are met. To simulate the effect of misregistration on QSM, we retrospectively performed an in-plane shift (10 mm in both x and y directions) on a middle slice from one stack of images that were acquired with consistent breath holds, and then computed the QSM.
Fat in the body produces artificial phase shifts, or chemical shift. If these chemical shifts are not removed prior to reconstructing QSM, the resulting QSM will have severe streaking artifacts. To accurately remove fat chemical shift, we first employed a magnitude image guide graph cuts method, for phase unwrapping and initial estimation of chemical shift. The joint field map estimation and water fat separation problem can be formulated as the following nonlinear least squares fitting problem of the complex MRI data with water and fat as the tissue chemical composition.
where s(TEn) is the measured MR signal at a certain echo time TEn, R2*=1/T2 is the decay rate, ƒF is the presumed fat chemical shift (ƒF≈−220 HZ for 1.5T or −440 HZ for 3T scanner), b is the field map (or field inhomogeneity), ρW is the water component, and ρF is the fat component. The above objective function is a non-convex function having multiple local minima. The solution to this minimization is especially sensitive to the initial guess of b. The complex signal generated by the two species can be formulated as
is the magnitude of s(TEn) and
is the phase of (ρw+ρFei2πƒ
{tilde over (ƒ)}=argminƒΣn=1N∥a(TEn)e−i2πbTE
For this nonlinear least square fitting problem, it is difficult to compute the analytical solutions, so here we provide an approximate analytical solution. Briefly, we use Taylor expansion to the first order near b and set the derivative of the residual function to zero to get the approximate solution
where
k is an integer number and k/ΔTE wrap is caused by the limited dynamic range. The true field map can be written as
b={tilde over (ƒ)}+k/ΔTE−mƒ
F−ψ(TEn,ρW,ρF), [46]
where m∈{0,1}, and the model error ψ(TEn, ρW, ρF) in the approximate solution is
Assuming the actual field map ƒs is spatially continuous and the model error is small compared to phase wraps or chemical shifts, appropriate k and m needs to be found to address the discontinuities in {tilde over (ƒ)}
k*,m*=argmink,m∥∇{tilde over (ƒ)}s∥22=argmink,m∥∇({tilde over (ƒ)}+k/ΔTE−mƒF)∥22 [48]
where ∇ denotes the gradient operator. This joint optimization over k and m can be efficiently solved using graph cuts with conditional jump moves. The model error in field map {tilde over (ƒ)}s is bounded by ψ(TEn, ρW, ρF). From the approximate solution of ψ(TEn, ρW, ρF), we find a greater last echo TEN leads to a smaller fitting error. In addition, ψ(TEn, ρW, ρF)=0 if ρW=0 or ρF=0. That is, if the voxels consist of pure water or pure fat, the model error of the fitting is zero. The true model error introduced in the single species fitting will be calculated in numerical simulation detailed in the methods section. It is verified that the true model error also decreases with increasing TEN and the model error ψ(TEn, ρW, ρF) is less than 15 Hz when TEN is larger than 17.6 ms, and the average error size is 2.67 Hz. Therefore, the above simultaneous phase unwrapping and chemical shift result b can be used as a good initial guess for IDEAL.
The above result was used in the above Eq.43 signal fitting to obtain water map and fat map. This simultaneous phase unwrapping and chemical shift determination based on graph cuts was used to initialize an iterative chemical shift update method to remove chemical shift from fat.
Since the T2 of lung is very short, there are usually no reliable signals from lung in GRE images. Thus, the myocardium of the heart is right next to region with no reliable signals, i.e. region with sources that generated background field. Current state of the art methods for removing or compensating for background field generally produce error in edges through either over estimation of background field, or erosion. To ensure an accurate suppression of the effect from background field in the final QSM, we used the Total Field Inversion method:
Where χ is the susceptibility map in the whole FOV, b is the total field generated by susceptibility estimated from the above fat-water decomposition, including susceptibility sources from tissue of interest and from background, d is the dipole kernel, ∇ is a gradient operator, λ is the regularization parameter, MG is the edge mask, and w is a SNR weighting. P is a preconditioner that increases the convergence speed of the algorithm. The preconditioner is a binary mask with 1 inside ROI and a larger value outside ROI. The outside value of P was empirically chosen to be 30.
Three of the five volunteers were able to perform adequate consistent breath holds for successful 2DBH QSM. In comparison, 3DNAV was successfully acquired on all five volunteers, and image quality of cardiac QSM from 3DNAV was consistently better than that from 2DBH on every subject.
As shown in
The average susceptibility difference between right and left ventricular blood pools, measured from the 2DBH QSM on the three volunteers who were able to perform adequate consistent breath holds, was 270±18 ppb, which corresponded to an oxygenation level of 79.6±1.1%. The average susceptibility difference between right and left ventricular blood pools, measured from the 3DNAV QSM on the five volunteers, was 262±13 ppb, which was 80±1% SvO2.
Our preliminary data showed the impracticality of the 2DBH approach for cardiac QSM. which was limited by its vulnerability to inconsistent breath holds. To understand this vulnerability, we simulated a single slice misregistration by shifting one of the slice in the volume that was successfully acquired by 2DBH approach, and we found drastically increased artifacts in QSM output. Therefore, as few as a single slice misregistration can cause major artifacts in the resulting QSM. The 3DNAV approach eliminated the need of breath holds, and the resulting cardiac QSM was uniformly better than cardiac QSM from 2DBH, making the 3DNAV approach a more viable option in clinical environment.
In summary, 3D acquisition with free-breathing navigator gating is a promising cardiac QSM approach without breath hold inconsistency and efforts.
In some cases, one or more of the above quantitative susceptibility mapping techniques can be implemented using the process 800 shown in
In some implementations, the process 800 can be used to transform MR signal data corresponding to a subject into susceptibility-based images that quantitatively depict the structure and/or composition of the subject. As an example, in some cases, the process 800 can be used to obtain MR data corresponding a subject, and process this MR data to generate a quantitative susceptibility map of the subject. Using this quantitative susceptibility map, one or more susceptibility-based images of the subject can be generated and displayed to a user. The user can then use these images for diagnostic or experimental purposes, for example to investigate the structure and/or composition and/or function of the subject, and/or to diagnose various conditions or diseases based, and/or to treat various conditions or diseases based, at least in part, on the images. As the process 800 can result in susceptibility-based images having higher quality and/or accuracy compared to other susceptibility mapping techniques, implementations of the process 800 can be used to improve a user's understanding of a subject's structure and/or composition and/or function, and can be used to improve the accuracy of any resulting medical diagnoses or experimental analyses.
The process 800 begins by acquiring magnetic resonance (MR) data corresponding to a subject (step 810). In some cases, the MR data can correspond to a patient (e.g., the entirety of a patient or a portion of a patient, such as a particular portion to the patient's body). In some cases, the MR data can correspond to one or more samples (e.g., an imaging phantom, a sample of one or more materials, a sample of one or more types of tissue, and/or one or more other objects).
In some implementations, the MR data can be acquired using an MRI scanner using one or more suitable pulse sequences. As an example, in some cases, MR data can be acquired using a gradient echo sequence that acquires MR data at a single echo time or at multiple echo times (e.g., two, three, four, five, and so forth). Various scan parameters can be used. As another example, MR data can be obtained using a 3D multi-echo spoiled gradient echo sequence on a 3T scanner (e.g., a GE Excite HD MR scanner) using a pulse sequence that samples at different echo times (TE) (e.g., 4 TEs, such as 3.8, 4.3, 6.3 and 16.3 ms) in an interleaved manner, and with following imaging parameters: repetition time (TR)=22 ms; voxel size=. 5×0.5×0.5 mm3, matrix size=256×64×64; bandwidth (BW)=±31.25 kHz; flip angle=15°. As another example, MR data can be obtained using a 3D multi-gradient echo sequence on a 3T scanner using imaging parameters TE/ΔTE/#TE/TR/FA/BW/FOV/matrix size=ms/5 ms/8/50 ms/20°/±62.50 kHz/21.8×24×12.8 cm3/232×256×64. Although an example sequences and example parameters are described above, these are merely illustrative. In practice, other sequences and parameters can be used, depending on various factors (e.g., the size of region to be examined, a known or assumed range of susceptibility values of the subject, scan time considerations, device limitations, and so forth).
After acquiring MR data corresponding to a subject, the process 800 continues by determining a magnetic field based on the MR data (step 820). As noted above, in some cases, the MRI signal phase is affected by the magnetic field corresponding to the magnetic susceptibility distribution of the subject and the chemical shift of tissue components.
The magnetic field can be determined from the complex MR data by fitting the detected signal as a sum over tissue spectral components, where each component signal characterized by an exponent with a negative real part representing signal decay and an imaginary part representing phase dependent on the magnetic field and the chemical shift (step 830). This fitting for such a complex signal model can be performed on iteratively using numerical optimization. A robust initialization for numerical optimization can be estimated using a simplified single species model and graph cut separating smooth field against chemical shift and phase unwrapping.
After determining the magnetic field based on the MR data, the process 800 continues by determining a relationship between the magnetic field and magnetic susceptibility (step 840). As described above, in some implementations, the relationship between the magnetic field and magnetic susceptibility can be expressed as a relationship between the magnetic field at a given location to the magnetic susceptibility at that location. In some cases, this relationship can be expressed in integral form, or in differential form. In the differential form, the relationship between the magnetic field at a given location to the magnetic susceptibility at that location can include an equation where a Laplacian of the magnetic field equals one third of the Laplacian of the magnetic susceptibility minus a second order derivative of the magnetic susceptibility.
After determining a relationship between the magnetic field and magnetic susceptibility, the process 800 continues by determining prior knowledge about tissue magnetic susceptibility distribution (step 850). One estimation of tissue susceptibility distribution is to use R2* values derived from the magnitude signal. High R2* values can be used to identify high susceptibility regions such as hemorrhages for preconditioning to accelerate the convergence of numerical optimization. An example of preconditioning is to separate the whole image volume into region with high susceptibility (including hemorrhages, air regions and background) and region with normal tissue susceptibility. Such preconditioning reduces the susceptibility search range and hence accelerates convergence. Another estimation of tissue susceptibility distribution is to use low (near zero) R2* regions for identifying regions of pure water such as cerebrospinal fluid in the ventricles in the brain, oxygenated arterial blood in the aorta, or regions of pure fat in the abdominal wall. As water and fat have known susceptibility values, they can be used to serve zero references (to water) to generate absolute susceptibility values using a minimal variance regularization.
After obtaining the susceptibility prior information and signal data noise property, the process 800 continues by estimating a magnetic susceptibility distribution of the subject based, at least in part, on the prior information and data noise property (step 860). As described above, estimating the susceptibility distribution of the subject can include determining a cost function corresponding to a susceptibility distribution, the magnetic field, masks corresponding to regions of interest. The cost function includes a data fidelity term based on data noise property expressing the relationship between the magnetic field and tissue susceptibility in an integral or differential form, and a regularization term expressing prior information. The estimated susceptibility distribution of the subject can be determined by identifying a particular susceptibility distribution that minimizes one or more of these cost functions. As described above, in some cases, this can be determined numerically.
After estimating a magnetic susceptibility distribution of the subject, the process 800 continues by generating one or more images of the subject based on the estimated susceptibility distribution of the subject (step 870). These images can be electronically displayed on a suitable display device (e.g., an LCD display device, LED display device, CRT display device, or other display device that can be configured to show images) and/or physically displayed on a suitable medium (e.g., printed, etched, painted, imprinted, or otherwise physically presented on a sheet or paper, plastic, or other material).
In some cases, the estimated susceptibility distribution of the subject can be visualized using a color scale, where each of several colors is mapped to a particular susceptibility value or range of susceptibility values. Accordingly, a two dimensional or three dimension image can be generated, where each pixel or voxel of the image corresponds to a particular spatial location of the subject, and the color of that pixel of voxel depicts the susceptibility value of the subject at that location.
In some cases, the color scale can include a gradient of colors and/or a gradient of gray shades (e.g., a gray scale), in order to depict a range of susceptibility values. In some cases, for a color scale that includes a gradient of colors or a gradient of gray shades, one end of the gradient can correspond to the lowest susceptibility value in a particular window of susceptibility values (e.g., an arbitrarily selected window of values), and the other end of the gradient can correspond to the highest susceptibility value in the window of values. For instance, for gray scale that includes a gradient of gray shades between pure white and pure black, pure white can be used to indicate the highest susceptibility value in a particular arbitrary window of values, and pure black can be used to indicate the lowest susceptibility value in the window of values. Other relationships between a color/gray scale and susceptibility values is also possible. For example, in some cases, for gray scale that includes a gradient of gray shades between pure white and pure black, pure white can be used to indicate the lower susceptibility value in a particular arbitrary window of values, and pure black can be used to indicate the highest susceptibility value in the window of values. Although the examples are described in the context of gray scales, in practice, similar relationships between color scales and susceptibility values are also possible.
Susceptibility values and colors can be mapped linearly (e.g., such that each absolute change in susceptibility value corresponds to a proportional change in color), logarithmically (e.g., such that each an exponential changes in susceptibility values correspond to a linear change in color), or according to any other mapping. Although examples mapping between color scales and susceptibility values are described above, these are merely illustrative examples. In practice, other mappings are also possible, depending on the implementation.
In some cases, the estimated susceptibility distribution can be regularized based on prior information regarding the subject. As an example, the estimated susceptibility distribution can regularized using one or more of the techniques described above. In some cases, regularizing the estimated susceptibility distribution based on prior information regarding the subject comprises determining boundary conditions associated with the subject.
Implementations of the above described techniques can be performed using a computer system.
The system 900 includes a processor 910, a memory 920, a storage device 930, and an input/output device 940. Each of the components 910, 920, 930, and 940 can be interconnected, for example, using a system bus 950. The processor 910 is capable of processing instructions for execution within the system 900. In some implementations, the processor 910 is a single-threaded processor. In some implementations, the processor 910 is a multi-threaded processor. In some implementations, the processor 910 is a quantum computer. The processor 910 is capable of processing instructions stored in the memory 920 or on the storage device 930. The processor 910 may execute operations such as performing one or more of the techniques described above.
The memory 920 stores information within the system 900. In some implementations, the memory 920 is a computer-readable medium. In some implementations, the memory 920 is a volatile memory unit. In some implementations, the memory 920 is a non-volatile memory unit.
The storage device 930 is capable of providing mass storage for the system 900. In some implementations, the storage device 930 is a non-transitory computer-readable medium. In various different implementations, the storage device 930 can include, for example, a hard disk device, an optical disk device, a solid-state drive, a flash drive, magnetic tape, or some other large capacity storage device. In some implementations, the storage device 930 may be a cloud storage device, e.g., a logical storage device including multiple physical storage devices distributed on a network and accessed using a network. In some examples, the storage device may store long-term data. The input/output device 940 provides input/output operations for the system 900. In some implementations, the input/output device 940 can include one or more of a network interface devices, e.g., an Ethernet card, a serial communication device, e.g., an RS-232 port, and/or a wireless interface device, e.g., an 802.11 card, a 3G wireless modem, a 4G wireless modem, etc. A network interface device allows the system 900 to communicate, for example, transmit and receive data. In some implementations, the input/output device can include driver devices configured to receive input data and send output data to other input/output devices, e.g., a keyboard, a mouse, a printer, a sensor (e.g., a sensor that measures component or system-related properties, a sensor that measures environmental-related properties, or other types of sensors), and a display device 960. In some implementations, mobile computing devices, mobile communication devices, and other devices can be used.
A computing system can be realized by instructions that upon execution cause one or more processing devices to carry out the processes and functions described above, for example, storing, maintaining, and displaying artifacts. Such instructions can include, for example, interpreted instructions such as script instructions, or executable code, or other instructions stored in a computer readable medium. A computing system can be distributively implemented over a network, such as a server farm, or a set of widely distributed servers or can be implemented in a single virtual device that includes multiple distributed devices that operate in coordination with one another. For example, one of the devices can control the other devices, or the devices may operate under a set of coordinated rules or protocols, or the devices may be coordinated in another fashion. The coordinated operation of the multiple distributed devices presents the appearance of operating as a single device.
Although an example processing system has been described in
The term “processing module” may encompass all apparatus, devices, and machines for processing data, including by way of example a programmable processor, a computer, or multiple processors or computers. A processing module can include, in addition to hardware, code that creates an execution environment for the computer program in question, e.g., code that constitutes processor firmware, a protocol stack, a database management system, an operating system, or a combination of one or more of them.
A computer program (also known as a program, software, software application, script, executable logic, or code) can be written in any form of programming language, including compiled or interpreted languages, or declarative or procedural languages, and it can be deployed in any form, including as a standalone program or as a module, component, subroutine, or other unit suitable for use in a computing environment. A computer program does not necessarily correspond to a file in a file system. A program can be stored in a portion of a file that holds other programs or data (e.g., one or more scripts stored in a markup language document), in a single file dedicated to the program in question, or in multiple coordinated files (e.g., files that store one or more modules, sub programs, or portions of code). A computer program can be deployed to be executed on one computer or on multiple computers that are located at one site or distributed across multiple sites and interconnected by a communication network.
Computer readable media suitable for storing computer program instructions and data include all forms of non-volatile or volatile memory, media and memory devices, including by way of example semiconductor memory devices, e.g., EPROM, EEPROM, and flash memory devices; magnetic disks, e.g., internal hard disks or removable disks or magnetic tapes; magneto optical disks; and CD-ROM and DVD-ROM disks. The processor and the memory can be supplemented by, or incorporated in, special purpose logic circuitry. Sometimes a server is a general purpose computer, and sometimes it is a custom-tailored special purpose electronic device, and sometimes it is a combination of these things. Implementations can include a back end component, e.g., a data server, or a middleware component, e.g., an application server, or a front end component, e.g., a client computer having a graphical user interface or a Web browser through which a user can interact with an implementation of the subject matter described is this specification, or any combination of one or more such back end, middleware, or front end components. The components of the system can be interconnected by any form or medium of digital data communication, e.g., a communication network. Examples of communication networks include a local area network (“LAN”) and a wide area network (“WAN”), e.g., the Internet.
Certain features that are described above in the context of separate implementations can also be implemented in combination in a single implementation. Conversely, features that are described in the context of a single implementation can be implemented in multiple implementations separately or in any sub-combinations.
The order in which operations are performed as described above can be altered. In certain circumstances, multitasking and parallel processing may be advantageous. The separation of system components in the implementations described above should not be understood as requiring such separation.
A number of embodiments of the invention have been described. Nevertheless, it will be understood that various modifications may be made without departing from the spirit and scope of the disclosure. Accordingly, other embodiments are within the scope of the following claims.
This application is a continuation application of U.S. application Ser. No. 17/142,475 filed on Jan. 6, 2021, which is a continuation application of U.S. application Ser. No. 15/943,753 filed on Apr. 3, 2018, now U.S. Pat. No. 10,890,641 issued on Jan. 12, 2021, which claims priority from U.S. Provisional Application No. 62/482,864, filed Apr. 7, 2017, the entire contents of which are incorporated herein by reference.
This invention was made with government support under Grant No. NS072370, NS090464, NS095562 and CA181566 awarded by the National Institutes of Health. The government has certain rights in this invention.
| Number | Date | Country | |
|---|---|---|---|
| 62482864 | Apr 2017 | US |
| Number | Date | Country | |
|---|---|---|---|
| Parent | 17142475 | Jan 2021 | US |
| Child | 18151860 | US | |
| Parent | 15943753 | Apr 2018 | US |
| Child | 17142475 | US |
