A General Framework for Automated Accurate Calculation of <i>b</i>-Matrix (Auto-<i>b</i>) in Diffusion MRI Pulse Sequences

Yuan, Lisha; Wu, Dan; He, Hongjian; Zhong, Jianhui

doi:https://doi.org/10.1155/2023/4610812

Concepts in Magnetic Resonance Part B, Magnetic Resonance Engineering

On this page

Abstract Introduction Discussion Conclusion Data Availability Conflicts of Interest Acknowledgments Supplementary Materials References Copyright Related Articles

Research Article | Open Access

Volume 2023 | Article ID 4610812 | https://doi.org/10.1155/2023/4610812

A General Framework for Automated Accurate Calculation of b-Matrix (Auto-b) in Diffusion MRI Pulse Sequences

Lisha Yuan,¹Dan Wu,¹Hongjian He,¹and Jianhui Zhong^1,2

Academic Editor: Lars G. Hanson

Received22 Aug 2022

Revised02 Mar 2023

Accepted17 Mar 2023

Published05 Apr 2023

Abstract

To derive accurate diffusion metrics, both imaging and diffusion-sensitizing gradient pulses should be accounted for when calculating the diffusion-weighted b-matrix. However, it is complex to derive analytical solutions due to complicated interactions between gradient pulses, including orthogonal directions. This study proposes a general framework to calculate the b-matrix automatically (dubbed as Auto-b). Based on the divide-and-conquer approach, the b-matrix calculation is appropriately segmented, and the symbolic mathematical library is applied to handle integration operations for each interval. If the specifications of all gradient pulses are provided to Auto-b, an accurate b-matrix can be obtained. Three examples are explored to validate the accuracy of Auto-b and to detect b-value errors when using approximate calculations. (1) In the conventional spin-echo example, Auto-b exhibits high accuracy, as indicated by the maximum relative deviation of 1.68‰ between its calculated b-matrices and those obtained from analytical expressions. (2) Auto-b is applied to investigate the contribution of imaging gradients to the b-matrix in an optimized spin-echo echo planar imaging sequence at submillimeter resolution. Specifically, ignoring the contribution of imaging gradients results in a b-value error of 12.16 s/mm² at the 0.8 × 0.8 × 0.8 mm³ resolution and 22.47 s/mm² at the 0.6 × 0.6 × 0.8 mm³ resolution, respectively, when nominal b = 0. (3) Auto-b is also utilized to analyze the influence of approximate calculations in the spatiotemporally encoded sequence. The results showed that neglecting the contribution of phase-encoding blips causes large b-value errors up to 11.02 s/mm². In addition, the rectangularization of trapezoidal waveforms led to a high relative b-value error of 39.91%. This study validates the high accuracy of Auto-b and underscores the importance of accurate b-value calculations in both submillimeter imaging and spatiotemporally encoded sequences. Attributed to its automation, accuracy, and broad applicability, Auto-b is helpful for developers of diffusion sequences.

1. Introduction

The diffusion process of water molecules in vivo can be detected to depict the characteristics and abnormalities of the microstructure of tissues [1, 2]. As a powerful noninvasive and quantitative technique, diffusion imaging has potential clinical utility in disease diagnosis, prognosis, and treatment evaluation in various human organs, such as the brain [3–6], spine [7–9], liver [10–12], breast [13, 14], and prostate [15, 16]. The b-matrix is a critical parameter in diffusion imaging as it summarizes the diffusion-induced signal attenuation determined by the direction, amplitude, and timing of all gradient pulses present in the sequence [17–19]. By analyzing the b-matrix and the detected signal, diffusion metrics can be derived.

When considering the effect of imaging gradients on diffusion-induced signal attenuation, conventional analytical approaches for calculating the b-matrix can become very complicated. In addition to diffusion-sensitizing gradients, various imaging gradients such as readout (RO), phase-encoding (PE), slice-selective, crusher, and possibly other gradients must be considered. The diversity and variety of gradient pulses make their interactions complicated and challenging to derive analytical solutions. On the one hand, each gradient needs to be split into small intervals [20–22], resulting in many intervals and a tedious deviation process. On the other hand, the analytical solution is only suitable for a specific sequence with specific arrangements of gradient pulses. Minor changes, such as a change of a gradient pulse in the shape or position, could invalidate the symbolic expression [20, 23].

One standard method to simplify the b-matrix calculation is reducing the contribution of imaging gradients to the b-matrix. Compared to conventional spin-echo (SE) [18, 23, 24] or echo planar imaging (EPI)-based [20] diffusion sequences, state-of-the-art diffusion sequences minimize the contribution of imaging gradients by placing the read-dephasing gradient close to RO gradients. In addition, crusher gradients, typically employed to eliminate interference from the free induction decay (FID) signal, are removed in the presence of diffusion-sensitizing gradients, as diffusion-sensitizing gradients can suppress the FID signal instead. As a result, the contribution of imaging gradients is generally ignored in the b-matrix calculation. However, at the submillimeter resolution, the effect of imaging gradients becomes a recurring issue. A high-performance gradient system facilitates in vivo human brain imaging with a submillimeter spatial resolution on 3T MR scanners [25–30]. Submillimeter diffusion MRI is gaining popularity because of its ability to characterize fine structural details, such as probing the short association fibers (U-fibers) [28]. As the amplitude of imaging gradients increases proportionally with increasing spatial resolution [18], it is unclear whether the contribution of imaging gradients remains negligible at submillimeter resolution.

An alternative method to reduce computational complexity is to apply approximate calculations, such as ignoring the contribution of RO gradients and PE blips [20] or treating trapezoidal waveforms as a rectangle [21, 24]. However, it is debatable whether errors introduced by approximate calculations are acceptable. While the EPI sequence suggests that PE blips have a minimal effect on diffusion-induced signal attenuation, caution should be taken with spatiotemporally encoded (SPEN) sequences due to the increased amplitude-duration product of PE blips [31]. In addition, it is essential to investigate errors associated with the approximate calculation based on the rectangularization of trapezoidal waveforms because of the quadratic dependence of the b-matrix on the gradient shape.

Considering the computational complexity and accuracy issues, this study proposes an automated framework to calculate the b-matrix (dubbed as Auto-b). To mitigate errors from approximate numerical integration and facilitate error analysis, the divide-and-conquer approach combined with the symbolic mathematical library has been applied to handle the integration operations in the b-matrix calculation. Given the specifications of all gradient pulses, our framework can convert each gradient into a series of subfunctions, break the whole time into proper intervals, calculate the b-matrix individually for each interval, and then add them to get the final result. A toolkit based on the divide-and-conquer approach is provided. Through comparison with analytical expressions, the accuracy of Auto-b was first validated with a conventional SE sequence. Auto-b was also applied to a currently optimized SE-EPI sequence to assess the contribution of imaging gradients to the b-value at submillimeter resolution. Finally, Auto-b was utilized in a SPEN sequence to investigate b-value errors introduced by approximate calculations.

2. Theory

2.1. A Generalized Definition of b-Matrix

Regarding sequences with multiple 180° refocusing pulses, the effective gradient [32] is adapted to capture the altered phase effect after rotations that change the sign of the precession. Given the original gradient time series , the effective gradient can be expressed as follows:where denotes the number of 180° refocusing pulses preceding time point and , , and are the components of the original gradient in the read, phase, and slice directions, respectively. It also applies to the case of sequences containing either zero or a single 180° refocusing pulse. Figure 1 illustrates the relationship between the original and effective gradient in the presence of five 180° refocusing pulses.

The b-matrix can then be expressed in terms of the effective gradient [22, 23] as follows:whereis an integral function summarizing the effect (i.e., area) of effective gradients along each axis up to time , and TE is the echo time. According to the vector multiplication rule, the element of the b-matrix can be obtained as follows:which demonstrates the interaction between gradient pulses along the same or orthogonal directions, and the indices i and j denote the coordinate direction of gradient pulses applied (i.e., the read, phase, or slice direction).

2.2. Key Technical Considerations for Automated b-Matrix Calculation

This study aims to generate the b-matrix result once the specifications of all the gradient pulses are given. To achieve this, the algorithm must handle the variations of gradient pulses across axes and sequences.

2.2.1. The Idea of Divide and Conquer

Due to the variations of the gradient pulses, it is challenging to obtain the analytical b-matrix expression directly for the whole time. To overcome this challenge, the idea of divide and conquer is adopted. We can break it down into smaller subproblems and solve each subproblem independently.

If we divide the whole period into N intervals, equation (4) can be adapted as follows:

Here, = 0 and = TE. The calculation of the b-matrix is then divided into two steps: (i) calculating in each interval separately; (ii) adding them to obtain the final value. Therefore, the main challenge is dividing the intervals appropriately and calculating the value in each interval.

2.2.2. Guidelines for Segmentation

Appropriate segmentation is essential to achieve divide and conquer. The whole time is generally segmented according to the gradient shape and antiphase instants (i.e., the center of the 180° refocusing pulse). Based on a schematic diagram of a pulse sequence (Figure 2), two guidelines for segmentation are described as follows:

(1) The Expression of Gradient Pulses in Each Interval Is a Single Expression. According to the waveform type, each gradient can be broken into intervals and expressed as a series of subfunctions of time. In practical terms, each interval must have only one algebraic expression. Taking two common waveform types (trapezoidal/sinusoidal) as an example, we describe how each gradient is segmented. Specifically, the first sinusoidal gradient along the x-axis can be expressed as

Meanwhile, the trapezoidal gradient along the y-axis can be expressed as three subfunctions:

Similarly, other waveforms that are separable into intervals can also yield time points for segmentation and be expressed as a series of subfunctions.

(2) The Subfunction Expression Should Not Cross Any Antiphase Instant. Since the expressions of effective gradients have an opposite sign before and after a 180° refocusing pulse, antiphase instants should also be included for segmentation. For example, the plateau period of the trapezoidal gradient along the y-axis should be divided into two segments:

This way the subsequent conversion from original gradient expressions to effective gradient expressions becomes convenient.

Following the two guidelines, we collect all the time points when either the expression of the gradient pulse or the expression meaning changes. According to the time points collected, the subfunctions in each interval can be constructed (see Table 1 as an example).

2.2.3. Calculating the b-Matrix at Each Interval

Once the whole time has been divided into proper intervals, the focus is on calculating the of each interval.

For any interval , the integral function can be expressed as

As long as is integrable within the interval, equation (9) becomeswhere is the indefinite integral of over the interval and . Thus, the symbolic mathematical library could be used to calculate . The same applies to the calculation of in this interval.

Similarly, if the product is integrable, then integrating the product yields the value of the interval (see equation (5)).

2.3. A General Framework for Automated and Accurate b-Matrix Calculation

This article proposes a general framework for implementing an automated b-matrix calculation (Auto-b). By providing gradient specifications as input (the first step), the subsequent steps of Auto-b can automatically calculate the b-matrix. Based on the divide-and-conquer approach, an accurate b-matrix can be obtained if a complete description of all gradient pulses is provided. Figure 3 depicts the flowchart of Auto-b, which consists of the following steps:

(1) Recording Gradient Specifications and Control Variables (User Input Required). The gradient specification includes its start time, shape (such as sinusoid or trapezoid), amplitude, ramp time, duration, repetition information, and direction. Specifically, the term “duration” represents the duration time of a gradient, except for a trapezoidal gradient shape, for which it represents the time from the initial rise to the end of its plateau. Table 2 shows an example of how gradient specifications can be documented (for instance, in an Excel spreadsheet) for gradient pulses shown in Figure 2.

Meanwhile, control variables consist of the excitation instant, the refocusing instant, and the antiphase instants. Note that the input of antiphase instants, which can be empty, a single value, or a group of values, varies depending on the number of 180° refocusing pulses in a sequence.

(2) Parsing Gradient Specifications. A specification item will be expanded if it contains duplicate gradients and all specification items on the same axis are combined. For instance, Table 3 lists all gradient pulses along the x-axis from Table 2. Each gradient is then resolved as a series of subfunctions, following the first guideline described in Section 2.2.2. Once all gradient pulses have been explained, a piecewise function of original gradients is formed for each axis.

(3) Constructing the Three-Dimensional (3D) Piecewise Function of the Original Gradients. To calculate the nondiagonal elements of the b-matrix, it is necessary to combine the piecewise function of the three axes. As mentioned in Section 2.2.2, the whole time is segmented based on the gradient shape and antiphase instants. Control variables are also utilized to restrict the time ranging from the excitation to the refocusing instant. By redefining the subfunctions of original gradients at each interval for three axes, a 3D piecewise function of original gradients is constructed (e.g., Table 1).

(4) Generating the Sequence Diagram. Based on the 3D piecewise function of original gradients, the sequence diagram is plotted and can be used to verify that gradient specifications and control variables are interpreted correctly.

(5) Calculating the b-Matrix Automatically. As the 3D piecewise function mentioned above is obtained for original gradients, the 3D piecewise function of effective gradients must first be obtained. Specifically, the subfunction of effective gradients for each interval is equal to or opposite to that of original gradients (equation (1)). Then, the value in each interval can be calculated, as described in Section 2.2.3, using the symbolic mathematical library. Finally, the values of all intervals are added to obtain the final b-matrix.

The Auto-b toolkit was developed based on the flowchart mentioned above and is accessible at https://github.com/lishayuan/calculate_b_matrix. The following is a brief description of the toolkit.

Firstly, the toolkit is capable of processing sequences with zero, single, or multiple 180° refocusing pulses. Demonstration examples of gradient-echo (GE), SE, SE-EPI, and turbo spin-echo (TSE) sequences are included to assist users in providing the correct input. Based on the toolkit, sequences with a single 180° refocusing pulse were evaluated in this study.

Secondly, the toolkit using the divide-and-conquer approach can theoretically handle any gradient shape that can be described analytically. Currently, only trapezoidal and sinusoidal waveforms have been implemented in the toolkit and evaluated in the article. As the toolkit is open source under the MIT license, users can easily incorporate a new gradient shape by defining its specification and adding the necessary code to convert it into subfunctions, similar to the trapezoidal and sinusoidal examples provided.

Finally, the toolkit has been enhanced to support numerical integration on a discrete representation of the waveforms, allowing the calculation of the b-matrix for more general waveforms. Please refer to the Discussion section for more details.

3. Tests of the Proposed Framework

In this study, three diffusion MRI sequences were utilized to investigate the accuracy of Auto-b, the contribution of imaging gradients at the submillimeter resolution, and the influence of approximate calculations on the b-value.

3.1. Accuracy of Auto-b for the Conventional SE Sequence

As described by Mattiello et al. [23], imaging gradients have a nonnegligible contribution to the b-matrix in the conventional SE sequences when the read-dephasing gradient is not placed close to the acquisition period or if there is a pair of crusher gradients (Figure 4(a)). A simulation protocol (Figure 4(b)) of the conventional SE sequence was used to verify the accuracy of Auto-b.

(a)

(b)

Figure 4

Illustration of a conventional SE imaging sequence (a) and a simulation protocol (b). In the simulation protocol, the variables i, t_i, G_i, and δ_i denote the type, start time, amplitude, and duration of the ith gradient pulse, respectively. The trapezoidal ramps have a uniform time of 200 μs. G_rdp represents the read-dephasing gradient; G_dr, G_dp, and G_ds represent diffusion-sensitizing gradients in the read, phase, and slice directions, respectively; G_cr, G_cp, and G_cs represent crusher gradients in the read, phase, and slice directions, respectively.

Table 4 exemplifies all the diagonal elements of the b-matrix calculated using Auto-b. As shown in Table 4, the value was 5.95 s/mm² when the crusher and diffusion-sensitizing gradients were absent, and the diagonal elements of the b-matrix increased with the increasing amplitude of crusher gradients. These confirm that the read-dephasing gradient and crusher gradients are essential in determining the b-matrix in the conventional SE sequence. Compared to the diagonal elements of b-matrices obtained from analytical expressions [23], the corresponding values calculated by Auto-b exhibited a maximum relative deviation of 1.68‰.

3.2. Contribution of Imaging Gradients at Submillimeter Resolution in an Optimized SE-EPI Sequence

In order to minimize the contribution of imaging gradients to the b-matrix, the SE-EPI sequence is optimized by placing the read-dephasing gradient close to the acquisition period (Figure 5) and avoiding applying crusher gradients when there is a pair of diffusion-sensitizing gradients (Figure 5(b)). However, a substantial contribution may still exist at submillimeter resolution. Drawing on imaging parameters employed for U-fiber imaging at 0.8 mm isotropic resolution [28], the following imaging protocols were used to investigate the contribution of imaging gradients at two submillimeter resolutions in the optimized SE-EPI sequence:

(a)

(b)

(1) Case 1 (0.8 × 0.8 × 0.8 mm³): matrix size = (RO) 256 × (PE) 112, TE = 58 ms, bandwidth = 1086 Hz/Px, and echo spacing = 1.01 ms; (2) Case 2 (0.6 × 0.6 × 0.8 mm³): matrix size = (RO) 342 × (PE) 150, TE = 81 ms, bandwidth = 860 Hz/Px, and echo spacing = 1.27 ms. Other imaging parameters were common, including a field of view (FOV) = (RO) 206 × (PE) 90 mm², TR = 3000 ms, 1 slice, partial Fourier in PE direction = 5/8, no in-plane acceleration, and with diffusion-weighting obtained along 3 orthogonal directions using 14 nominal b-values (b = 0, 10, 20, 30, 40, 50, 100, 200, 300, 400, 500, 600, 700, and 800 s/mm²). The sequence timing was generated using a 3T MR scanner with a maximum gradient amplitude of 80 mT/m and a maximum slew rate of 200 T/m/s. The nominal b-value in this study only considers the contribution of diffusion-sensitizing gradients, while the accurate b-value calculated by Auto-b also includes the contribution of imaging gradients.

The discrepancy between accurate and nominal b-values was first examined (Figure 6(a)). For nominal b = 0, imaging gradients contributed 12.16 s/mm² and 22.47 s/mm² to the b-value for Cases 1 and 2, respectively. When the nominal b-value changes from 100 s/mm² to 800 s/mm², ignoring the contribution of imaging gradients led to b-value errors ranging from 1.04 s/mm² to 5.38 s/mm² for Case 1 and from 1.20 s/mm² to 4.74 s/mm² for Case 2. Moreover, the effect of b-value errors on the apparent diffusion coefficient (ADC) was investigated by calculating the ADC value using the data acquired at nominal b = 0 and the other b-value (Figure 6(b)). As the nominal b-value was selected from 800 s/mm² to 100 s/mm², the relative error of ADC increased from 0.85% to 9.96% for Case 1 and from 2.22% to 20.32% for Case 2.

(a)

(b)

Considering the accuracy of small b-values may profoundly affect the evaluation of physiological parameters in the intravoxel incoherent motion (IVIM) model, a simulation experiment was conducted. The IVIM model, which enables simultaneous evaluation of tissue diffusion and perfusion, can be described by the biexponential formula [33]:where S represents the signal intensity in the presence of diffusion-sensitizing gradients, S₀ represents the signal intensity in the absence of diffusion-sensitizing gradients, f represents the perfusion fraction, and D and D^∗ are the diffusion and pseudodiffusion coefficients, respectively. Assuming ideal physiological parameters (S₀ = 1000, D = 2.2 × 10⁻³ mm²/s, f = 0.3, and D^∗ = 1.6 × 10⁻² mm²/s), signals were simulated with accurate b-values based on equation (11). Physiological parameters were then estimated using a 2-step fitting procedure [34, 35]. Table 5 shows that fitting with accurate b-values resulted in a negligible error for all physiological parameters. On the other hand, fitting with nominal b-values resulted in a relative error of 29.05% for the pseudodiffusion coefficient at the 0.8 × 0.8 × 0.8 mm³ resolution and as high as 46.91% at the 0.6 × 0.6 × 0.8 mm³ resolution.

3.3. Influence of Approximate Calculations in the SPEN Sequence

SPEN methods have received much attention due to their ability to greatly reduce distortions in geometry and intensity [36, 37]. To briefly illustrate the principle, Figure 7(a) shows an axial sketch of human kidneys to be acquired (with the y-axis as the SPEN axis), and Figure 7(b) shows a SPEN sequence. During the 90° chirp pulse, positions are excited sequentially along the SPEN axis. Given the chirp pulse with a symmetric linearly decreasing frequency and a simultaneous encoding gradient with a positive polarity, positions at the positive edge of the SPEN axis are excited first, and positions at the negative edge are excited last. During acquisition, positions along the SPEN axis were refocused in reverse order, resulting in spatially dependent diffusion weighting [21, 38]. Here is the definition of the b-matrix for the SPEN sequence,where and represent the unique excitation and refocusing instants for each position .

(a)

(b)

To investigate the impact of approximate calculations on the b-matrix in the SPEN sequence, two commonly used approximations were assessed: (1) the approximation based on a uniform refocusing instant and (2) the rectangularization of trapezoidal waveforms. Specifically, the investigation focused on gradient pulses along the SPEN axis, and both accurate and approximate values were calculated for the analysis.

3.3.1. The Approximation Based on a Uniform Refocusing Instant

The b-matrix calculation for the SPEN sequence can be accurately performed by assigning a specific refocusing instant to different positions [31, 36], as shown in Figure 7(b). However, the approximation based on a uniform refocusing instant assumes that all positions are refocused at the center instant of the last acquired echo, which ignores diffusion-induced signal attenuation from PE blips. To investigate the impact of this approximation, the following imaging protocols were set up with different FOV scenarios:

(1) Scenario 1: FOV = (RO) 220 × (PE) 220 mm², matrix size = 128 × 128, TE = 118 ms, bandwidth = 1502 Hz/Px, and echo spacing = 0.75 ms; (2) Scenario 2: FOV = (RO) 220 × (PE) 110 mm², matrix size = 128 × 64, TE = 84 ms, bandwidth = 1302 Hz/Px, and echo spacing = 0.93 ms; (3) Scenario 3: FOV = (RO) 220 × (PE) 66 mm², matrix size = 220 × 64, TE = 85 ms, bandwidth = 1336 Hz/Px, and echo spacing = 0.93 ms. The three scenarios shared common image parameters, including TR = 3000 ms, FA = 90°, 5.0 mm slice thickness, 1 slice, with diffusion-weighting obtained along the SPEN axis using 5 nominal b-values (b = 0, 200, 400, 600, and 800 s/mm²), and a chirp pulse with a bandwidth of 50 kHz and a duration of 5.12 ms.

Figure 8(a) illustrates three imaging FOVs selected on the human kidney sketch, and Figure 8(b) displays the corresponding errors along the SPEN axis for each FOV scenario. Since the contribution of PE blips to the b-value is independent of diffusion-sensitizing gradients, the result of Figure 8(b) is consistent for different nominal b-values. However, the impact of the approximation based on a uniform refocusing instant varies depending on the position along the SPEN axis. Positions at the positive edge were not unaffected, while positions at the negative edge exhibited the largest error. Specifically, as the FOV along the SPEN axis changed from 220 mm and 110 mm to 66 mm, the absolute error for positions at the negative edge increased from 1.67 s/mm² and 4.05 s/mm², to 11.02 s/mm², respectively.

(a)

(b)

3.3.2. Rectangularization of Trapezoidal Waveforms

The imaging protocol of Scenario 3 was further utilized to investigate errors associated with the rectangularization of trapezoidal waveforms. To simplify the evaluation, only diffusion-sensitizing gradients were rectangularized. Three possible rectangularization schemes were discussed:

(1) Rect 1, which discarded the ramps but kept the amplitude of diffusion-sensitizing gradients (Figure 9(b)); (2) Rect 2, which kept the base and area with the amplitude adjusted (Figure 9(c)); and (3) Rect 3, which kept the area and amplitude with the base adjusted (Figure 9(d)).

Table 6 represents the values and their relative error associated with the rectangularization of diffusion-sensitizing gradients at both edges of the SPEN axis. For positions at the negative edge, the accurate values were comparable to nominal values. In contrast, approximate calculations based on Rect 1, Rect 2, and Rect 3 exhibited a relative error of 34.29%, 1.29%, and −0.13%, respectively, for all nominal b-values. For position at the positive edge, due to the opposite polarity between the encoding and diffusion-sensitizing gradients, the accurate values were lower than the nominal values. In addition, the relative errors produced by each rectangularization scheme were similar but slightly greater than those observed at the negative edge position.

4. Discussion

Auto-b has been proposed for an automated and accurate calculation of the b-matrix. The accuracy of Auto-b was validated in the conventional SE sequence by comparing b-matrices calculated by Auto-b with those obtained from analytical expressions. Moreover, the SE-EPI and SPEN examples demonstrated the importance of accurate b-matrix calculation in sequences with high-resolution or strong imaging gradients.

As a general framework, Auto-b can be applied to various sequences and gradient shapes. Although only sequences with a single 180° refocusing pulse are evaluated in this study, the b-matrix in terms of the effective gradient can be calculated for sequences containing zero, single, or multiple 180° refocusing pulses. This allows Auto-b to cover a wide range of sequences, such as GE, SE, double SE, SE-EPI, and TSE sequences. Based on the divide-and-conquer approach, Auto-b theoretically supports a wide range of waveform types, including trapezoidal and sinusoidal waveforms, as well as other waveforms that can be described by polynomial, trigonometric, exponential, or rational functions. In addition, Auto-b has been enhanced to support numerical integration, allowing the b-matrix to be calculated for more general waveforms. After discretizing the waveform or providing discrete gradient data directly (especially for waveforms that cannot be described analytically), the b-matrix can be calculated using the numerical integration approach (see equations S1-S2 in the Supplementary Materials for numerical integration on a discrete representation of the waveforms).

The SE example demonstrated the high accuracy of Auto-b, with a maximum relative deviation of 1.68‰ between its calculated b-matrices and those obtained from analytical expressions. Auto-b calculates the b-matrix using a divide-and-conquer approach combined with a symbolic mathematical library, ideally producing identical results to analytical expressions. The observed deviation suggests that a large cumulative rounding error has been found in the results obtained from the implemented toolkit. The potential source of the rounding error could be equation (10), which calculates and stores and as numerical values for the integration operation at each interval. The error accumulates over the whole time and eventually affects the accuracy. In contrast, the analytical solution involves numerical calculations only when the values are substituted into the final analytical expression, resulting in a minor rounding error.

Despite optimizing the SE-EPI sequence, neglecting the contribution of imaging gradients to the b-matrix causes large errors in the b-value and estimated diffusion metrics. A notable contribution of imaging gradients was first found in submillimeter imaging. Specifically, crusher gradients substantially contribute to the b-value at nominal b = 0, followed by the interaction between the slice-selective and diffusion-sensitizing gradients at nominal b ≠ 0. Since the moment of the crusher gradients in this study is determined by the voxel size in the read direction, their contribution is proportional to the in-plane resolution. In addition, the interaction between the slice-selective and diffusion-sensitizing gradients increases with increasing nominal b-value or decreasing voxel size in the slice direction. If crusher gradients are present at nominal b ≠ 0 and interact with the diffusion-sensitizing gradients, the accurate b-value will further deviate from the nominal one. The SE-EPI example also demonstrates the importance of a precise b-matrix in estimating diffusion metrics. Due to the large b-value deviation at nominal b = 0, the relative error of the ADC is greater when using the other data acquired at a small b-value. Assuming that the pseudodiffusion coefficient is much larger than the diffusion coefficient, the pseudodiffusion effect can only be detected at small b-values [34, 35]. In the IVIM model, the estimation of the pseudodiffusion coefficient is greatly affected by the high relative error of small b-values. In contrast, the diffusion coefficient, determined by high b-values (greater than 300 s/mm² in this study), produces negligible errors at both submillimeter resolutions.

Some limitations need to be acknowledged. As the separation of data and algorithm eliminates the need for in-depth knowledge of the b-matrix definition or programming language, it is user-friendly to use the toolkit to perform the calculation. However, the information about gradient specifications may only be available to some customers, and time is required for users to prepare the input data. A higher level of automation can be achieved if vendors provide a convenient interface for customers to access this information. In addition, Auto-b currently supports only the case of a single coherence pathway involving multiple 180° refocusing pulses, and sequences refocusing multiple coherence pathways (such as TSE variants) are excluded. The b-matrix definition in terms of the effective gradient applies only to refocusing pulses with a flip angle of 180°, and sequences containing refocusing pulses with other flip angles, implying multiple coherence pathways, are beyond the scope of this article. Finally, in addition to imaging gradients discussed in this study, other confounding factors such as gradient nonlinearity [39, 40], concomitant fields [41], eddy currents [42], and gradient miscalibration [43, 44] may also cause the actual b-matrix to deviate from the nominal value. Although the alteration of the gradient waveform due to the imperfect gradient field is reduced by the hardware upgrade and improved sequence design, a comprehensive evaluation is still required in future work.

5. Conclusion

Auto-b is proposed to calculate the b-matrix for diffusion sequences automatically. Once all the specifications of gradient pulses are provided to Auto-b, an accurate b-matrix can be obtained automatically. Auto-b, based on the divide-and-conquer approach, has been shown to reproduce the analytical results of a conventional SE sequence accurately and highlights the importance of accurate b-matrix calculation for sequences with high resolution or strong imaging gradients. Due to its automation and accuracy, Auto-b helps developers to calculate the b-matrix for various diffusion sequences. In future studies, in vivo experiments will be designed to further investigate the validity of Auto-b in accurately estimating diffusion metrics.

Data Availability

The toolkit implemented in MATLAB is available at https://github.com/lishayuan/calculate_b_matrix. It supports both the divide-and-conquer and numerical integration approaches and includes demonstration examples of GE, SE, SE-EPI, TSE, and SPEN sequences.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

The authors would like to thank Congbo Cai and Dehe Weng for helpful suggestions on the experimental design and Lixia Yuan and Wei Liu for careful comments on this article. This work was supported by the Fundamental Research Funds of the Central Universities (Grant no. 2021FZZX002-19).

Supplementary Materials

Auto-b is extended to support the approach based on numerical integration, which enables the b-matrix to be calculated for more general waveforms. (Supplementary Materials)

References

D. Le Bihan and M. Iima, “Diffusion magnetic resonance imaging: what water tells us about biological tissues,” PLoS Biology, vol. 13, no. 7, Article ID 1002203, 2015.
View at: Publisher Site | Google Scholar
P. B. Kingsley, “Introduction to diffusion tensor imaging mathematics: Part I. Tensors, rotations, and eigenvectors,” Concepts in Magnetic Resonance: Part A, vol. 28, no. 2, pp. 101–122, 2006.
View at: Publisher Site | Google Scholar
Y. Assaf and O. Pasternak, “Diffusion tensor imaging (DTI)-based white matter mapping in brain research: a review,” Journal of Molecular Neuroscience, vol. 34, no. 1, pp. 51–61, 2008.
View at: Publisher Site | Google Scholar
J. Goveas, L. O’Dwyer, M. Mascalchi et al., “Diffusion-MRI in neurodegenerative disorders,” Magnetic Resonance Imaging, vol. 33, no. 7, pp. 853–876, 2015.
View at: Publisher Site | Google Scholar
A. Delouche, A. Attyé, O. Heck et al., “Diffusion MRI: pitfalls, literature review and future directions of research in mild traumatic brain injury,” European Journal of Radiology, vol. 85, no. 1, pp. 25–30, 2016.
View at: Publisher Site | Google Scholar
M. Ouyang, J. Dubois, Q. Yu, P. Mukherjee, and H. Huang, “Delineation of early brain development from fetuses to infants with diffusion MRI and beyond,” NeuroImage, vol. 185, pp. 836–850, 2019.
View at: Publisher Site | Google Scholar
M. Cercignani, M. A. Horsfield, F. Agosta, and M. Filippi, “Sensitivity-encoded diffusion tensor MR imaging of the cervical cord,” AJNR. American journal of neuroradiology, vol. 24, no. 6, pp. 1254–1256, 2003.
View at: Google Scholar
Y. Eguchi, S. Ohtori, M. Yamashita et al., “Diffusion magnetic resonance imaging to differentiate degenerative from infectious endplate abnormalities in the lumbar spine,” Spine, vol. 36, no. 3, pp. E198–E202, 2011.
View at: Publisher Site | Google Scholar
O. A. Zaninovich, M. J. Avila, M. Kay, J. L. Becker, R. J. Hurlbert, and N. L. Martirosyan, “The role of diffusion tensor imaging in the diagnosis, prognosis, and assessment of recovery and treatment of spinal cord injury: a systematic review,” Neurosurgical Focus, vol. 46, no. 3, p. E7, 2019.
View at: Publisher Site | Google Scholar
B. Taouli, A. J. Tolia, M. Losada et al., “Diffusion-weighted MRI for quantification of liver fibrosis: preliminary experience,” American Journal of Roentgenology, vol. 189, no. 4, pp. 799–806, 2007.
View at: Publisher Site | Google Scholar
J. Zhang, Y. Guo, X. Tan et al., “MRI-based estimation of liver function by intravoxel incoherent motion diffusion-weighted imaging,” Magnetic Resonance Imaging, vol. 34, no. 8, pp. 1220–1225, 2016.
View at: Publisher Site | Google Scholar
M. M. McNamara, J. V. Thomas, L. F. Alexander et al., “Diffusion-weighted MRI as a screening tool for hepatocellular carcinoma in cirrhotic livers: correlation with explant data—a pilot study,” Abdominal Radiology, vol. 43, no. 10, pp. 2686–2692, 2018.
View at: Publisher Site | Google Scholar
E. Yılmaz, O. Sarı, A. Yılmaz et al., “Diffusion-weighted imaging for the discrimination of benign and malignant breast masses; utility of ADC and relative ADC,” Journal of the Belgian Society of Radiology, vol. 102, no. 1, pp. 24–26, 2018.
View at: Publisher Site | Google Scholar
N. Amornsiripanitch, S. Bickelhaupt, H. J. Shin et al., “Diffusion-weighted MRI for unenhanced breast cancer screening,” Radiology, vol. 293, no. 3, pp. 504–520, 2019.
View at: Publisher Site | Google Scholar
G. Lemberskiy, A. B. Rosenkrantz, J. Veraart, S. S. Taneja, D. S. Novikov, and E. Fieremans, “Time-dependent diffusion in prostate cancer,” Investigative Radiology, vol. 52, no. 7, pp. 405–411, 2017.
View at: Publisher Site | Google Scholar
X. Yang, L. Shan, H. Cao, X. Jiang, and X. Ma, “A valuable MRI examination method for prostate cancer screening,” Medicine, vol. 99, no. 48, Article ID e23134, 2020.
View at: Publisher Site | Google Scholar
P. J. Basser, J. Mattiello, and D. Lebihan, “Estimation of the effective self-diffusion tensor from the NMR spin echo,” Journal of Magnetic Resonance - Series B, vol. 103, no. 3, pp. 247–254, 1994.
View at: Publisher Site | Google Scholar
D. Güllmar, J. Haueisen, and J. R. Reichenbach, “Analysis of b-value calculations in diffusion weighted and diffusion tensor imaging,” Concepts in Magnetic Resonance: Part A, vol. 25, no. 1, pp. 53–66, 2005.
View at: Publisher Site | Google Scholar
P. B. Kingsley, “Introduction to diffusion tensor imaging mathematics: Part II. Anisotropy, diffusion‐weighting factors, and gradient encoding schemes,” Concepts in Magnetic Resonance: Part A, vol. 28, no. 2, pp. 123–154, 2006.
View at: Publisher Site | Google Scholar
J. Mattiello, P. J. Basser, and D. Le Bihan, “The b matrix in diffusion tensor echo-planar imaging,” Magnetic Resonance in Medicine, vol. 37, no. 2, pp. 292–300, 1997.
View at: Publisher Site | Google Scholar
E. Solomon, N. Shemesh, and L. Frydman, “Diffusion weighted MRI by spatiotemporal encoding: analytical description and in vivo validations,” Journal of Magnetic Resonance, vol. 232, pp. 76–86, 2013.
View at: Publisher Site | Google Scholar
M. Zubkov, T. Stait-Gardner, and W. S. Price, “Efficient and precise calculation of the b-matrix elements in diffusion-weighted imaging pulse sequences,” Journal of Magnetic Resonance, vol. 243, pp. 65–73, 2014.
View at: Publisher Site | Google Scholar
J. Mattiello, P. J. Basser, and D. LeBihan, “Analytical expressions for the b matrix in NMR diffusion imaging and spectroscopy,” Journal of Magnetic Resonance - Series A, vol. 108, no. 2, pp. 131–141, 1994.
View at: Publisher Site | Google Scholar
O. Brihuega-Moreno, F. P. Heese, C. Tejos, and L. D. Hall, “Effects of, and corrections for, cross‐term interactions in q-space MRI,” Magnetic Resonance in Medicine, vol. 51, no. 5, pp. 1048–1054, 2004.
View at: Publisher Site | Google Scholar
H. C. Chang, M. Sundman, L. Petit et al., “Human brain diffusion tensor imaging at submillimeter isotropic resolution on a 3 Tesla clinical MRI scanner,” NeuroImage, vol. 118, pp. 667–675, 2015.
View at: Publisher Site | Google Scholar
K. Setsompop, Q. Fan, J. Stockmann et al., “High‐resolution in vivo diffusion imaging of the human brain with generalized slice dithered enhanced resolution: simultaneous multislice (gSlider‐SMS),” Magnetic Resonance in Medicine, vol. 79, no. 1, pp. 141–151, 2018.
View at: Publisher Site | Google Scholar
G. Ramos-Llordén, L. Ning, C. Liao et al., “High-fidelity, accelerated whole-brain submillimeter in vivo diffusion MRI using gSlider-spherical ridgelets (gSlider-SR),” Magnetic Resonance in Medicine, vol. 84, no. 4, pp. 1781–1795, 2020.
View at: Publisher Site | Google Scholar
F. Movahedian Attar, E. Kirilina, D. Haenelt et al., “Mapping short association fibers in the early cortical visual processing stream using in vivo diffusion tractography,” Cerebral Cortex, vol. 30, no. 8, pp. 4496–4514, 2020.
View at: Publisher Site | Google Scholar
Y. Hu, X. Wang, Q. Tian et al., “Multi-shotdiffusion-weighted MRI reconstruction with magnitude-basedspatial-angular locally low-rank regularization (SPA-LLR),” Magnetic Resonance in Medicine, vol. 83, no. 5, pp. 1596–1607, 2020.
View at: Publisher Site | Google Scholar
G. Ramos-Llordén, G. Vegas-Sánchez-Ferrero, C. Liao, C. Westin, K. Setsompop, and Y. Rathi, “SNR-enhanced diffusion MRI with structure-preservinglow-rank denoising in reproducing kernel Hilbert spaces,” Magnetic Resonance in Medicine, vol. 86, no. 3, pp. 1614–1632, 2021.
View at: Publisher Site | Google Scholar
Y. Chen, J. Li, X. Qu et al., “Partial Fourier transform reconstruction for single-shot MRI with linear frequency-swept excitation,” Magnetic Resonance in Medicine, vol. 69, no. 5, pp. 1326–1336, 2013.
View at: Publisher Site | Google Scholar
R. Karlicek and I. J. Lowe, “A modified pulsed gradient technique for measuring diffusion in the presence of large background gradients,” Journal of Magnetic Resonance, vol. 37, no. 1, pp. 75–91, 1980.
View at: Publisher Site | Google Scholar
D. Le Bihan, E. Breton, D. Lallemand, M. L. Aubin, J. Vignaud, and M. Laval-Jeantet, “Separation of diffusion and perfusion in intravoxel incoherent motion MR imaging,” Radiology, vol. 168, no. 2, pp. 497–505, 1988.
View at: Publisher Site | Google Scholar
C. Federau, P. Maeder, K. O’Brien, P. Browaeys, R. Meuli, and P. Hagmann, “Quantitative measurement of brain perfusion with intravoxel incoherent motion MR imaging,” Radiology, vol. 265, no. 3, pp. 874–881, 2012.
View at: Publisher Site | Google Scholar
C. Federau, R. Meuli, K. O’Brien, P. Maeder, and P. Hagmann, “Perfusion measurement in brain gliomas with intravoxel incoherent motion MRI,” American Journal of Neuroradiology, vol. 35, no. 2, pp. 256–262, 2014.
View at: Publisher Site | Google Scholar
A. Tal and L. Frydman, “Spatial encoding and the single-scan acquisition of high definition MR images in inhomogeneous fields,” Journal of Magnetic Resonance, vol. 182, no. 2, pp. 179–194, 2006.
View at: Publisher Site | Google Scholar
N. Ben-Eliezer, Y. Shrot, and L. Frydman, “High-definition, single-scan 2D MRI in inhomogeneous fields using spatial encoding methods,” Magnetic Resonance Imaging, vol. 28, no. 1, pp. 77–86, 2010.
View at: Publisher Site | Google Scholar
E. Solomon, G. Liberman, N. Nissan, and L. Frydman, “Robust diffusion tensor imaging by spatiotemporal encoding: principles and in vivo demonstrations,” Magnetic Resonance in Medicine, vol. 77, no. 3, pp. 1124–1133, 2016.
View at: Publisher Site | Google Scholar
H. Y. Mesri, S. David, M. A. Viergever, and A. Leemans, “The adverse effect of gradient nonlinearities on diffusion MRI: from voxels to group studies,” NeuroImage, vol. 205, Article ID 116127, 2020.
View at: Publisher Site | Google Scholar
D. I. Malyarenko, B. D. Ross, and T. L. Chenevert, “Analysis and correction of gradient nonlinearity bias in apparent diffusion coefficient measurements,” Magnetic Resonance in Medicine, vol. 71, no. 3, pp. 1312–1323, 2014.
View at: Publisher Site | Google Scholar
M. A. Bernstein, X. J. Zhou, J. A. Polzin et al., “Concomitant gradient terms in phase contrast MR: analysis and correction,” Magnetic Resonance in Medicine, vol. 39, no. 2, pp. 300–308, 1998.
View at: Publisher Site | Google Scholar
R. W. Chan, C. von Deuster, D. Giese et al., “Characterization and correction of eddy-current artifacts in unipolar and bipolar diffusion sequences using magnetic field monitoring,” Journal of Magnetic Resonance, vol. 244, pp. 74–84, 2014.
View at: Publisher Site | Google Scholar
S. Mohammadi, Z. Nagy, H. E. Möller et al., “The effect of local perturbation fields on human DTI: characterisation, measurement and correction,” NeuroImage, vol. 60, no. 1, pp. 562–570, 2012.
View at: Publisher Site | Google Scholar
Y. Lee, A. O. Kettinger, B. J. Wilm et al., “A comprehensive approach for correcting voxel‐wise b‐value errors in diffusion MRI,” Magnetic Resonance in Medicine, vol. 83, no. 6, pp. 2173–2184, 2020.
View at: Publisher Site | Google Scholar

Copyright

Copyright © 2023 Lisha Yuan et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

469

Downloads

295

Citations