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

,


Introduction
Te difusion 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, difusion imaging has potential clinical utility in disease diagnosis, prognosis, and treatment evaluation in various human organs, such as the brain [3][4][5][6], spine [7][8][9], liver [10][11][12], breast [13,14], and prostate [15,16]. Te b-matrix is a critical parameter in difusion imaging as it summarizes the difusion-induced signal attenuation determined by the direction, amplitude, and timing of all gradient pulses present in the sequence [17][18][19]. By analyzing the b-matrix and the detected signal, difusion metrics can be derived.
When considering the efect of imaging gradients on difusion-induced signal attenuation, conventional analytical approaches for calculating the b-matrix can become very complicated. In addition to difusion-sensitizing gradients, various imaging gradients such as readout (RO), phaseencoding (PE), slice-selective, crusher, and possibly other gradients must be considered. Te 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][21][22], resulting in many intervals and a tedious deviation process.
On the other hand, the analytical solution is only suitable for a specifc sequence with specifc 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 difusion 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 difusion-sensitizing gradients, as difusion-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 efect 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][26][27][28][29][30]. Submillimeter difusion MRI is gaining popularity because of its ability to characterize fne structural details, such as probing the short association fbers (U-fbers) [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 efect on difusion-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 specifcations 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 fnal result. A toolkit based on the divide-and-conquer approach is provided. Trough comparison with analytical expressions, the accuracy of Autob was frst validated with a conventional SE sequence. Autob 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.

A Generalized Defnition of b-Matrix.
Regarding sequences with multiple 180°refocusing pulses, the efective gradient [32] is adapted to capture the altered phase efect after rotations that change the sign of the precession. Given the original gradient time series , the efective gradient can be expressed as follows: where k denotes the number of 180°refocusing pulses preceding time point t and G r (t), G p (t), and G s (t) 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 efective gradient in the presence of fve 180°r efocusing pulses. Te b-matrix can then be expressed in terms of the efective gradient [22,23] as follows: where is an integral function summarizing the efect (i.e., area) of efective gradients along each axis up to time t, 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).

Key Technical Considerations for Automated b-Matrix
Calculation. Tis study aims to generate the b-matrix result once the specifcations of all the gradient pulses are given. To achieve this, the algorithm must handle the variations of gradient pulses across axes and sequences.

Te 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, t 0 � 0 and t N � TE. Te calculation of the b-matrix is then divided into two steps: (i) calculating b ij in each interval separately; (ii) adding them to obtain the fnal value. Terefore, the main challenge is dividing the intervals appropriately and calculating the b ij value in each interval.

Guidelines for Segmentation.
Appropriate segmentation is essential to achieve divide and conquer. Te 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) Te 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. Specifcally, the frst 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) Te Subfunction Expression Should Not Cross Any Antiphase Instant. Since the expressions of efective 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: 180º 90º   Tis way the subsequent conversion from original gradient expressions to efective 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).

Calculating the b-Matrix at Each
Interval. Once the whole time has been divided into proper intervals, the focus is on calculating the b ij of each interval. For As long as G i (t) is integrable within the interval, equation (9) becomes Tus, the symbolic mathematical library could be used to calculate F i (t). Te same applies to the calculation of F j (t) in this interval.
Similarly, if the product F i (t)•F j (t) is integrable, then integrating the product yields the b ij value of the interval (see equation (5)).

A General Framework for Automated and Accurate b-Matrix Calculation.
Tis article proposes a general framework for implementing an automated b-matrix calculation (Auto-b). By providing gradient specifcations as input (the frst 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 fowchart of Auto-b, which consists of the following steps: (1) Recording Gradient Specifcations and Control Variables (User Input Required). Te gradient specifcation includes its start time, shape (such as sinusoid or trapezoid), amplitude, ramp time, duration, repetition information, and direction. Specifcally, 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 specifcations 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 Specifcations. A specifcation item will be expanded if it contains duplicate gradients and all specifcation 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 frst 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 Tree-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 redefning 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 specifcations 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 efective gradients must frst be obtained. Specifcally, the subfunction of Table 1: Te piecewise function of gradient pulses in Figure 2.
Concepts in Magnetic Resonance Part B, Magnetic Resonance Engineering 5 efective gradients for each interval is equal to or opposite to that of original gradients (equation (1)). Ten, the b ij value in each interval can be calculated, as described in Section 2.2.3, using the symbolic mathematical library. Finally, the b ij values of all intervals are added to obtain the fnal b-matrix. Te Auto-b toolkit was developed based on the fowchart mentioned above and is accessible at https://github.com/ lishayuan/calculate_b_matrix. Te 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 defning its specifcation 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.

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

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. Table 4 exemplifes all the diagonal elements of the bmatrix calculated using Auto-b. As shown in Table 4, the b rr value was 5.95 s/mm 2 when the crusher and difusionsensitizing gradients were absent, and the diagonal elements of the b-matrix increased with the increasing amplitude of crusher gradients. Tese confrm that the readdephasing 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‰.

Contribution of Imaging Gradients at Submillimeter
Resolution in an Optimized SE-EPI Sequence. In order to minimize the contribution of imaging gradients to the bmatrix, 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 difusion-sensitizing gradients (Figure 5(b)). However, a substantial contribution may still exist at submillimeter resolution. Drawing on imaging parameters employed for U-fber 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: ( Te discrepancy between accurate and nominal b-values was frst examined (Figure 6(a)). For nominal b � 0, imaging gradients contributed 12.16 s/mm 2 and 22.47 s/mm 2 to the b-value for Cases 1 and 2, respectively. When the nominal bvalue changes from 100 s/mm 2 to 800 s/mm 2 , ignoring the contribution of imaging gradients led to b-value errors ranging from 1.04 s/mm 2 to 5.38 s/mm 2 for Case 1 and from 1.20 s/mm 2 to 4.74 s/mm 2 for Case 2. Moreover, the efect of b-value errors on the apparent difusion coefcient (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 2 to 100 s/mm 2 , 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.
Considering the accuracy of small b-values may profoundly afect the evaluation of physiological parameters in the intravoxel incoherent motion (IVIM) model, where S represents the signal intensity in the presence of difusion-sensitizing gradients, S 0 represents the signal intensity in the absence of difusion-sensitizing gradients, f represents the perfusion fraction, and D and D * are the difusion and pseudodifusion coefcients, respectively. Assuming ideal physiological parameters (S 0 = 1000, D = 2.2 × 10 − 3 mm 2 /s, f = 0.3, and D * = 1.6 × 10 − 2 mm 2 /s), signals were simulated with accurate b-values based on equation (11). Physiological parameters were then estimated using a 2-step ftting procedure [34,35]. Table 5 shows that ftting with accurate b-values resulted in a negligible error for all physiological parameters. On the other hand, ftting with nominal b-values resulted in a relative error of 29.05% for the pseudodifusion coefcient at the 0.8 × 0.8 × 0.8 mm 3 resolution and as high as 46.91% at the 0.6 × 0.6 × 0.8 mm 3 resolution.

Infuence 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 briefy 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 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. Te trapezoidal ramps have a uniform time of 200 μs. G rdp represents the read-dephasing gradient; G dr , G dp , and G ds represent difusion-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. where t e (y) and t a (y) represent the unique excitation and refocusing instants for each position y.
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. Specifcally, the investigation focused on gradient pulses along the SPEN axis, and both accurate and approximate b yy values were calculated for the analysis.

Te Approximation Based on a Uniform Refocusing
Instant. Te b-matrix calculation for the SPEN sequence can be accurately performed by assigning a specifc refocusing instant to diferent 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 difusion-induced signal attenuation from PE blips. To investigate the impact of this approximation, the following imaging protocols were set up with diferent FOV scenarios: (  Figure 8(a) illustrates three imaging FOVs selected on the human kidney sketch, and Figure 8(b) displays the corresponding b yy errors along the SPEN axis for each FOV scenario. Since the contribution of PE blips to the b-value is independent of difusion-sensitizing gradients, the result of Figure 8 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 unafected, while positions at the negative edge exhibited the largest error. Specifcally, as the FOV along the SPEN axis changed from 220 mm and 110 mm to 66 mm, the absolute b yy error for positions at the negative edge increased from 1.67 s/mm 2 and 4.05 s/mm 2 , to 11.02 s/mm 2 , respectively.

Rectangularization of Trapezoidal Waveforms.
Te imaging protocol of Scenario 3 was further utilized to investigate errors associated with the rectangularization of trapezoidal waveforms. To simplify the evaluation, only difusion-sensitizing gradients were rectangularized. Tree possible rectangularization schemes were discussed: (1) Rect 1, which discarded the ramps but kept the amplitude of difusion-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 b yy values and their relative error associated with the rectangularization of difusionsensitizing gradients at both edges of the SPEN axis. For positions at the negative edge, the accurate b yy 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 difusion-sensitizing gradients, the accurate b yy 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.

Discussion
Auto-b has been proposed for an automated and accurate calculation of the b-matrix. Te accuracy of Auto-b was     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 efective gradient can be calculated for sequences containing zero, single, or multiple 180°refocusing pulses. Tis 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 bmatrix 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).
Te SE example demonstrated the high accuracy of Autob, 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 divideand-conquer approach combined with a symbolic mathematical library, ideally producing identical results to analytical expressions. Te observed deviation suggests that a large cumulative rounding error has been found in the results obtained from the implemented toolkit. Te potential source of the rounding error could be equation (10), which calculates and stores F i (t s− 1 ) and h i (t s− 1 ) as numerical values for the integration operation at each interval. Te error accumulates over the whole time and eventually afects the accuracy. In contrast, the analytical solution involves numerical calculations only when the values are substituted into the fnal 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 difusion metrics. A notable contribution of imaging gradients was frst found in submillimeter imaging. Specifcally, crusher gradients substantially contribute to the b-value at nominal b = 0, followed by the interaction between the slice-selective and difusionsensitizing 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 difusion-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 difusionsensitizing gradients, the accurate b-value will further deviate from the nominal one. Te SE-EPI example also demonstrates the importance of a precise b-matrix in estimating difusion 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 pseudodifusion coefcient is much larger than the difusion coefcient, the pseudodifusion efect can only be detected at small b-values [34,35]. In the IVIM model, the estimation of the pseudodifusion coefcient is greatly afected by the high relative error of small b-values. In contrast, the difusion coefcient, determined by high bvalues (greater than 300 s/mm 2 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 defnition or programming language, it is user-friendly to use the toolkit to perform the calculation. However, the information about gradient specifcations 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. Te b-matrix defnition in terms of the efective gradient applies only to refocusing pulses with a fip angle of 180°, and sequences containing refocusing pulses with other fip 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 felds [41], eddy currents [42], and gradient miscalibration [43,44] may also cause the actual bmatrix to deviate from the nominal value. Although the alteration of the gradient waveform due to the imperfect gradient feld is reduced by the hardware upgrade and improved sequence design, a comprehensive evaluation is still required in future work.

Conclusion
Auto-b is proposed to calculate the b-matrix for difusion sequences automatically. Once all the specifcations of gradient pulses are provided to Auto-b, an accurate b-matrix can be obtained automatically. Auto-b, based on the divideand-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 difusion sequences. In future studies, in vivo experiments will be designed to further investigate the validity of Auto-b in accurately estimating difusion metrics.

Data Availability
Te toolkit implemented in MATLAB is available at https:// github.com/lishayuan/calculate_b_matrix. It supports both Concepts in Magnetic Resonance Part B, Magnetic Resonance Engineering the divide-and-conquer and numerical integration approaches and includes demonstration examples of GE, SE, SE-EPI, TSE, and SPEN sequences.