Investigation of Human Cardiac Fiber Architecture with Multiple diffusion MRI schemes at Multi-Scales based on Simulation Data

Background: Diffusion tensor imaging (DTI), diffusion spectrum imaging (DSI) and Q-ball imaging (QBI) are currently three main diffusion MRI (dMRI) schemes available for non-invasive investigation of cardiac fiber architecture. Although DSI and QBI have undoubtedly greater potential to reveal complex cardiac fiber structures than DTI, it however remains unclear to which level and at which scale they provide more gain for investigating cardiac fiber structure. Method: This work attends to provide a quantitative description of cardiac fiber architecture derived from different schemes at various scales. Due to the limit of the spatial resolution of clinical MRI scanner and with the absence of the ground-truth, it is difficult to give the accurate description. To deal with this issue, we simulate firstly DTI, DSI and QBI of a cardiac fiber model with the structure a priori known at different scales, and then the estimation accuracy, the diffusion metrics and the helix and transverse angles of cardiac fiber obtained by different schemes at different scales are calculated. Results: The results show that although DSI and QBI can distinguish multiple fiber orientations, they are readily to generate false positive and false negative fibers which influence therefore the estimation accuracy. When there are multiple fiber orientations in one voxel, the diffusion anisotropy detected by DTI is higher than DSI and QBI, the range of helix and transverse angle decreases with increasing of the scales, and that detected by DSI is larger than DTI and QBI. Conclusion: The results showed that the proposed dMRI simulator provides a valuable tool for simulating realistic DW images of whole human hearts, which can be used as the gold standard to study the fiber structures of the heart.

4 Fig. 2. DW images along one diffusion direction obtained with DTI, DSI and QBI at three scales. The corresponding b-values are 1000 mm 2 /s, 1843 mm 2 /s and 3000mm 2 /s.

Estimation accuracy
The estimation accuracy in number of cardiac fiber orientations for DTI, DSI and QBI at three scales was given in Table 2. It can be observed that the performance of all the schemes was dropping with the increase of the voxel size. At scale 1 where was present one single fiber orientation inside a voxel, DTI and DSI were able to estimate accurately the number of orientations, but the accuracy of QBI was only about 98%. At scale 2, the accuracy of DSI is higher but that of QBI was lower than DTI. At scale 3 where may exist multiple orientations inside a voxel the performance of both DSI and QBI is better than that of DTI. Moreover, at scales 2 and 3, the false negative error for DSI and QBI was higher than false positive, which means that the number of fiber orientations in most of voxels was underestimated.  Fractional anisotropy Fig. 4 shows the visual and quantitative description of diffusion anisotropies produced by different schemes at three scales.
It observed that the diffusion anisotropy decreased with the increasing of the scales for all the schemes. At scale 1, the mean diffusion anisotropy obtained by DSI was larger than that of DTI, which in turn was larger than QBI. However, the STD of the 6 diffusion anisotropy of QBI was the greatest at this scale, followed by DSI and then DTI. At scale 2, DSI and QBI had the similar diffusion anisotropy maps, however, DTI beard rather higher anisotropy values. At scale 3, the variation of diffusion anisotropy from epicardium to endocardium turned to be more obvious for DSI but that for DTI and QBI were not evident, and the diffusion anisotropy in the midcardium obtained by DTI was larger than DSI which was larger than QBI.

Helix and transverse angles
Fig . 5 shows the cardiac fiber orientations of one short-axis slice obtained by DTI, DSI and QBI at three scales. The fiber orientations are represented as false color-coded helix angle (HA) (Fig. 5(a)) and transverse angle (TA) maps ( Fig. 5(b)). These helix and transverse angles range from -90° to 90°, where the HA of -90° (or 90°) represents a fiber pointing toward the base (apex) of the heart, and the TA of -90º (or 90°) means a fiber pointing toward (outward) the ventricular cavity. We observed that, visually, the change of the imaging schemes did not influence greatly the measured HA and TA maps which were almost the same for DTI, DSI and QBI. To quantitatively compare the fiber orientations obtained by different dMRI schemes, a short-axis slice of LV was divided into six myocardial region of interests (ROIs) corresponding to American Heart Association (AHA) segments [51], and then the joint histograms between (transmural) location and (helix and transverse) angle for DTI, DSI and QBI at three scales in 6 AHA segments were calculated. In Fig. 6, each column of the joint histogram gives the distribution of angles at a given transmural location and each line the distribution of transmural locations at a given angle. The number in top left corner of the histogram is the correlation factor (CF) that indicates the linear dependence of the helix or transverse angle on the location.  We observed that, the CF of HA obtained by DSI and QBI was smaller than DTI at all the scales and it did not change so much with the variation of scales; CF of TA decreased in AHA 4 and AHA 5 but increased in AHA 1 and AHA 3 with the increase of imaging scales for all imaging schemes. The variation of CF of TA in AHA 2 and AHA 6 was not regular. We can also see that the probability of a given HA or TA value at a specified location varied along with the imaging schemes.
For AHA 1, the probability of larger HA values at mid-endocardium and endocardium obtained by QBI was greater than DSI and DTI at the three scales. Especially at scale 1, the probability of HA larger than 36° obtained by DSI and QBI was zero but that obtained by QBI was about 1%, which resulted in that QBI had a relative larger mean HA value at mid-endocardium and endocardium, as illustrated in Fig. 7(a). At epicardium, DSI and QBI had a greater probability to experience a HA between -18° and 0°. At midcardium, the difference between three schemes was not obvious at scale 1, however, at scale 2, the probability of HA around zero for DSI was lower than that of DTI and QBI. At scale 3, the HA values at epicardium generated by DSI were smaller than DTI and QBI.
For AHA 2, the difference between the three schemes at scale 1 was not obvious, however, at scale 2 and 3, such difference became significant. We observed that, DSI and QBI generated a higher probability to have smaller HA values at epicardium.
At mid-endocardium, the HA value obtained by DSI was lower than that of DTI and QBI.
For AHA 3, at scale 1, the joint histogram for different schemes were similar. At scale 2 and 3, the HA at midcardium for DSI had a greater probability to be positive than DTI and QBI.
For AHA 4, the joint histogram of DTI and DSI were similar but that of QBI was a little different from them, at scale 1 and scale 2, the main difference was located at mid-endocardium, where QBI has relatively greater HA values than DTI and DSI.
At scale 3, the difference at epicardium turned to be more obvious, where the HA value obtained by QBI was the lowest.
For AHA 5, the main difference between three schemes located at the mid-epicardium, where DSI had larger probability than DTI and QBI to generate bigger HA values.
For AHA 6, the distributions of joint histogram obtained by DSI at epicardium, midcardium and mid-endocardium was different from QBI and DTI. At scale 1, the range of HAs at epicardium obtained by DSI was narrower than DTI and QBI. At scale 2 and scale 3, the joint histogram at midcardium experienced an obvious change through the three schemes.
The joint histogram distribution of TA was illustrated in Fig. 6(b). TA and HA have similar variation trend. The main difference between the three schemes occurred at AHA 2 (epicardium at scale 2 and scale 3), AHA 4 and AHA 6 (epicardium and endocardium).
In order to describe quantitatively the variation of HA and TA values, the left ventricle was divided into five layers from epicardium to endocardium and the mean±std of HA and TA values for each layer of all the AHA segments were then calculated, as illustrated in Fig. 7. We can observe that, at the same imaging scale, the main difference among the three schemes in terms of mean±std value of HA and TA located at the zones of AHA 2 (epicardium), AHA 5 (mid-epicardium), and AHA 6 (midepicardium). In addition, such mean±std varies clearly with imaging scale. As regard to HA values, the range of which at AHA 10 1, AHA 2, and AHA 4 decreased with increasing of imaging scales and at the other zones increased.

Discussions
In this work, we investigated the cardiac fiber architecture using DTI, DSI and QBI at three scales, the quantitative description was given in terms of estimation accuracy, the diffusion anisotropy, and the helix and transverse angles.
From the point of view of estimation accuracy of cardiac fiber structure, we found that the performance of DSI and QBI at scale 3 was better than DTI as expected, but the improvement was not evident. Although they are able to distinguish multiple fiber, they are readily to overestimate or underestimate the cardiac fiber numbers due to the limit of the reconstruction methods.
At scale 1 and scale 2, the estimation accuracy of QBI was inferior to DTI. This can be explained by the nature of QBI that reconstructs the ODF (orientation distribution function) based on spherical harmonics (SH), in which the number of diffusion directions determines the maximum SH order that can be used in the reconstruction, and then truncating SH order will smooth the ODF. Such smoothing will influence the performance of QBI, by introducing spurious or fibers, decreasing angular resolution and increasing the deviation angle between the estimated and true fiber orientations, as illustrated in Fig. 3. This finding is also reported in the works of Wilkins [52]and Schilling [53]. In addition, we noticed that, with increase of the scale, the false positive error decreased and the false negative error increased, which implies that DSI and QBI may be able to distinguish two crossing orientations accurately but not several ones.
As to the diffusion anisotropy, since QBI generates many spurious fibers with small fractional anisotropy, which correspondingly leads to the minimum mean value of nQA. In contrast, the mean nQA obtained by DSI is larger than DTI at scale 1, which may be caused by the large b-value. In our DSI simulation, the maximum b-value is 4000 s/mm 2 that corresponds to a q value of 54.9 mm -1 , such large q value is able to produce a smaller detected fiber diameter, as reported in the work of Huang [54], and therefore a larger anisotropy indices. This finding is consistent with the reports of Schilling [53] and Mastropietro [55], in which they found that large b values lead to large diffusion anisotropy. In addition, the standard deviation of nQA obtained by DSI and QBI is much greater than that of FA, which verifies that the diffusion metrics of DSI and QBI can better reflect the heterogeneity of fiber orientations. This can also be observed in Fig. 4(a), the nQA of DSI and QBI has a ring form distribution at scale 2 and scale 3. At midcardium, the fiber orientation heterogeneity is smaller and therefore nQA is larger, but at endocardium or epicardium this is reversed.
In the light of cardiac fiber orientations reconstructed by different schemes, we observed that CF of TA values changed obviously with imaging scales at segments AHA 1, AHA 2, AHA 4 and AHA 5. Such finding may be caused by the fiber crossing. Generally, mean±std of HA and TA maps for different AHA segments are a common index to quantitatively describe cardiac fiber orientations [12,13]. However, it has not yet been reported that the mean±std of fiber angles varies with imaging scales and reconstruction methods. Our findings showed that the range of HA at AHA1, AHA 2 and AHA 3 and that of TA at AHA 1, AHA2, AHA4, and AHA 5 decreased dramatically with the increase of imaging scales. This demonstrates that we should take into account the effect of imaging voxel sizes when comparing the HA or TA ranges of cardiac fiber orientations.
Besides, we observed that the difference in mean±std of HA and TA maps between different imaging schemes were not obvious.
However, the variation of joint histogram between fiber angle and transmural location can indicate the difference between schemes in detail. In our work, the resolution of joint histogram was selected as 10x10. The bigger the window size of joint histogram, the more details it gives. Since the cardiac fiber crossing or the heterogeneity of fiber orientations in a voxel is not evident, the averaging process is easy to eliminating the difference between different schemes. This is why we can tell the difference between schemes from joint histogram rather than mean±std maps.
In short, due to the nature of cardiac fiber structure and that of DSI and QBI reconstruction algorithms, the performance of DSI and QBI is not always better than DTI. Their performance is greatly dependent on the cardiac fiber configurations (crossing angle), the b-values used, the order of SH basis and the noise level, etc.

Conclusion
We have investigated cardiac fiber structure at different imaging scales with DTI, DSI and QBI respectively. The schemes were compared in terms of the estimation of the number of fibers, orientations accuracy, diffusion metrics, helix and transvers angles of fibers. Our experimental results demonstrated that, DSI and QBI are able to distinguish multiple fiber orientations but are readily to generate the false positive and false negative errors. The diffusion anisotropy detected by DTI is higher than DSI and QBI when there are multiple fiber orientations in one voxel. At the zones where may present the crossing fibers, the range of helix and transverse angle decreases with increasing of the scales, and that detected by DSI is larger than DTI and QBI.

Diffusion MRI data simulation
To obtain the simulated DW images, firstly a virtual cardiac fiber model was created. The dataset used to model the cardiac fiber structure originated from the Center of cardiovascular bioinformatics and modeling of Johns Hopkins University [45]. It is under the form of vector fields with a spatial resolution of 0.43×0.43×1mm 3 and image size of 256×256×134. This dataset allows us to obtain the fiber orientation distribution of an entire human heart in 3D. For the cardiac fiber modeling, we assumed that, at each voxel, there were several cardiac myocytes modeled by oriented cylinders with a fix direction. The length of cylinders varies from 50 to 100 μm and their diameter ranges from 10 to 20 μm. Once the cardiac fiber structure was obtained, the Mont-Carlo method was used to simulate the diffusion process of water molecules in such structure and the DW images were calculated based on dMRI theory. The detailed modeling and simulation process can be referred to our previous work [46].
To save simulation time and computing resources, DW images for different diffusion encoding schemes at three scales were simulated with a spatial resolution of 0.43×0.43×1, 0.86×0.86×1 and 1.72×1.72×1 mm 3 respectively. The detailed simulation parameters were summarized in Table 1. The diffusion encoding directions of DTI and QBI were distributed uniformly over a sphere, but those of DSI were sampled from the Cartesian q-space with the minimum and maximum q-values of 15 and 55 mm -1 respectively.

Processing of diffusion MRI simulation data
For DTI data, second-order symmetric tensors were reconstructed at each voxel [20], from which the local cardiac fiber orientation and fractional anisotropy (FA) [18], were extracted: where 1  , 2  and 3  are the eigenvalues of tensor. The eigenvector corresponding to the maximum eigenvalue points to the local fiber orientation.
For DSI data, the reconstruction was achieved with the following steps [24]: the diffusion spectrum was first calculated by taking the 3D Fourier transform of the simulated q-space DW images. A Hanning filter with size of 17 was pre-multiplied to the diffusion spectrum images to smooth signal at high q value. Then, diffusion ODF deconvolution [47] was conducted using a regularization parameter of 0.1. From such ODF, the fiber orientations were finally extracted by searching local maximum.
For QBI data, the ODF was reconstructed using spherical harmonic-based method [48] with SH order of 8 and the diffusion ODF deconvolution was conducted using a regularization parameter of 0.006. In both DSI and QBI schemes, the diffusion index frequently used was the generalized fractional anisotropy (GFA) and normalized quantitative anisotropy (nQA) defined as [40]  iso  is isotropic background diffusion ODF, and 0 Z is a scaling constant that scales free water diffusion to 1 and makes QA value comparable across subjects. All the above reconstructions and fiber tracking were implemented with DSI studio. (http://dsi-studio.labsolver.org/ ).
Helix and transverse angles are two most important parameters to describe the cardiac fiber architecture, which were calculated using a local cylindrical coordinate system that is formed by l e , c e and r e [49,50] as depicted in Fig. 1.

Evaluation criteria
Two criteria were used to evaluate the accuracy of reconstructions of different schemes at various scales, one is related to the fiber number in each voxel, which expressed by the right ratio, false positive and false negative errors, and another is the reconstruction accuracy of fiber orientations. Right ratio is defined as the proportion of voxels in which a reconstruction algorithm estimates correctly the number of fiber orientations. False positive error means the estimated fiber number is more than the real one, and false negative is contrary to that. 16 The ground-truth for the fiber orientation and number of orientations in each voxel at different scales is given by the cardiac fiber model. During the modeling process, we assumed that at each voxel of size 0.43 × 0.43 × 1 mm 3 is present one cardiac fiber orientation, the ground-truth for the number of cardiac fiber orientation is therefore one at such scale. As the voxel size increased to 0.86 × 0.86 × 1 mm 3 and 1.72 × 1.72 × 1 mm 3 , the ground-truth for the number of orientations and the corresponding fiber orientations in each voxel is estimated as follows: when the angle between each two orientations is smaller than a threshold ( In this work is set to 6°), we assumed that they have the same orientation. Based on this assumption, the cardiac fibers are divided into several groups, the number of groups is defined therefore as the number of orientations in this voxel, and mean fiber orientation of each group are taken as the ground truth for fiber orientation.
The estimation accuracy in fiber orientation is assessed in terms of the mean deviation angle between the estimated orientation and the true one

Availability of data and material
The datasets analyzed during the current study available from the corresponding author on reasonable request.

Ethics approval and consent to participate
Not applicable.