Probing cardiomyocyte mobility with multi-phase cardiac diffusion tensor MRI

Purpose Cardiomyocyte organization and performance underlie cardiac function, but the in vivo mobility of these cells during contraction and filling remains difficult to probe. Herein, a novel trigger delay (TD) scout sequence was used to acquire high in-plane resolution (1.6 mm) Spin-Echo (SE) cardiac diffusion tensor imaging (cDTI) at three distinct cardiac phases. The objective was to characterize cardiomyocyte organization and mobility throughout the cardiac cycle in healthy volunteers. Materials and methods Nine healthy volunteers were imaged with cDTI at three distinct cardiac phases (early systole, late systole, and diastasis). The sequence used a free-breathing Spin-Echo (SE) cDTI protocol (b-values = 350s/mm2, twelve diffusion encoding directions, eight repetitions) to acquire high-resolution images (1.6x1.6x8mm3) at 3T in ~7 minutes/cardiac phase. Helix Angle (HA), Helix Angle Range (HAR), E2 angle (E2A), Transverse Angle (TA), Mean Diffusivity (MD), diffusion tensor eigenvalues (λ1-2-3), and Fractional Anisotropy (FA) in the left ventricle (LV) were characterized. Results Images from the patient-specific TD scout sequence demonstrated that SE cDTI acquisition was possible at early systole, late systole, and diastasis in 78%, 100% and 67% of the cases, respectively. At the mid-ventricular level, mobility (reported as median [IQR]) was observed in HAR between early systole and late systole (76.9 [72.6, 80.5]° vs 96.6 [85.9, 100.3]°, p<0.001). E2A also changed significantly between early systole, late systole, and diastasis (27.7 [20.8, 35.1]° vs 45.2 [42.1, 49]° vs 20.7 [16.6, 26.4]°, p<0.001). Conclusion We demonstrate that it is possible to probe cardiomyocyte mobility using multi-phase and high resolution cDTI. In healthy volunteers, aggregate cardiomyocytes re-orient themselves more longitudinally during contraction, while cardiomyocyte sheetlets tilt radially during wall thickening. These observations provide new insights into the three-dimensional mobility of myocardial microstructure during systolic contraction.

Introduction the sequence parameters and the gradient design [13,14]. The second limitation inherent to M1-M2 compensated gradient waveforms is that it is assumed that the underlying cardiac motion is comprised of only constant velocity or acceleration during encoding. Consequently, M1-M2 compensated gradients lead to intravoxel signal loss when cardiac motion is more complex. In general, these assumptions may only be true during a subset of cardiac phases, which differ from patient to patient or even beat to beat. This subject-specific dependence results in 'on-the-fly' manual adjustments to find the correct trigger delay (TD) for high quality imaging, adding to the overall acquisition time and increasing protocol complexity.
The objective of this work was to acquire SE high-resolution cDTI data at multiple cardiac phases to characterize cardiomyocyte organization and mobility at several cardiac phases. A prospective TD scout approach was used to estimate the patient specific cardiac phases available for cDTI. Two to three distinct cardiac phases were acquired in healthy volunteers in which in vivo cardiomyocyte mobility during contraction and filling was characterized.

Imaging protocol
Nine healthy volunteers (N = 9) were recruited for an MRI exam after obtaining signed statements of informed consent under an IRB and HIPPA compliant protocol approved by the University of California Los Angeles (UCLA) ethic committee (IRB#14-001743). All acquisitions were performed on a single 3T MRI scanner (Prisma, Siemens), using an 18-channel body phased-array coil and a 32-channel spine array coil. Long-axis, four-chamber, and shortaxis localizers, as well as short-axis balanced steady-state free precession (bSSFP) cine images were first acquired under breath-hold conditions. A short-axis cine DENSE sequence (resolution 2.3 x 2.3 mm, 3D four-point displacement encoding, k e = 0.8, TE = 1.05ms, TR = 15ms, 30 cardiac phases or more) was also performed using free-breathing and an expiratory triggered navigator for a subset of subjects (N = 5). The measured displacements were used to estimate the longitudinal cardiac displacement between the cardiac phases of interest.

cDTI sequence
The cDTI acquisition parameters were optimized to minimize signal dropout induced by cardiac motion and achieve higher spatial resolution. A custom single-shot spin-echo EPI sequence was modified to incorporate second-order (M 1 -M 2 ) motion compensated diffusion encoding gradients [14,15], inner-volume reduced-field-of-view (rFOV) [17], prospective navigator echo for respiratory triggering (acceptance windows 4 mm during expiration) and fat saturation. First, the sequence parameters were chosen so that the combination of the rFOV, the readout bandwidth, partial Fourier factor, and GRAPPA [18] parallel imaging minimized the EPI echo train length. The resultant short EPI echo train length limited the image distortion caused by high resolution and motion blurring. This approach provides a minimal TE contribution from the readout (Fig 1C), which permits a longer time for diffusion encoding, and reduces both the overall TE and diffusion encoding temporal footprint of the sequence. Second, the EPI readout distortion was further reduced by using a cardiac B 0 shim focused on the left ventricle. Lastly, the M 1 -M 2 gradient waveform was optimized on-the-fly symmetrically to reach a second gradient moment (M 2 ) as close to zero as possible as defined by Welsh et al. [13].
For each subject, one mid-ventricular short-axis slice was imaged at different cardiac phases. The acquisition parameters were: TE = 61 ms, rFOV 200x160 mm 2 , acquisition matrix 128x104, resolution 1.6x1.6x8 mm 3 (interpolated in-plane to 0.8x0.8x8 mm 3 ), parallel imaging with 2x-GRAPPA, partial Fourier 6/8, and readout bandwidth of 1860 Hz/px. Two b-values (0 and 350 s/mm 2 ) were acquired using twelve diffusion encoding directions and eight repetitions. The diffusion encoding duration with motion compensation was 48 ms while the EPI readout was 23 ms. After optimization, the four diffusion encoding gradients have total durations of 7.2/12.6/12.6/7.2 ms with a ramp time of 1.52 ms, corresponding to gradient moment The TD scout was acquired during free-breathing with cardiac synchronization. Images were acquired every two heartbeats to allow for signal recovery (TR�2s). (B) The cDTI acquisitions was cardiac and respiratory triggered. One image was acquired every respiratory cycle (TR�4s). (C) Details of the cDTI sequence used in this work. The timings of the EPI acquisition T nav /δ 90˚/ δ 180˚/ T a /T b /T EPI-echo /T EPI /T E were 8/2/4/24/20/9/23/61ms and the timings of the repeated diffusion encoding gradient δ 1 /δ 2 were 7.2/12.6ms. The deadtime δ d was 4ms.
https://doi.org/10.1371/journal.pone.0241996.g001 M 0 = 0 mT/m/s, M 1 = 1.3 mT/m/s 2 , M 2 = 62.6 mT/m/s 3 with an amplitude G = 71.3 mT/m (for reference a monopolar gradient waveform with the same sequence parameters has M 0 = 0 mT/m/s, M 1 = 11.5 x 10 3 mT/m/s 2 , M 2 = 555.9 x 10 3 mT/m/s 3 and G = 22.5 mT/m). The single-shot acquisition was ECG and end-respiratory triggered ( Fig 1B) using a navigator located on the right liver dome resulting in a TR of one breathing cycle (~4000 ms). A total of 104 images were obtained per cardiac phase corresponding to an acquisition time of~7 minutes.

TD scout acquisitions
To precisely determine which cardiac phases were suited for cardiac diffusion imaging, a prospective TD scout acquisition was implemented [15,19]. As shown in Fig 1A, the TD scout consists of repeated diffusion acquisitions where the TD between the electrocardiogram Rwave and the excitation pulse of the diffusion sequence is sequentially shifted at each acquisition. The TD scout was acquired during free-breathing without navigator triggering and one image was acquired every two heartbeats (TR � 2000ms for a heart rate of 60 bpm). A total of thirty cardiac phases were scouted over the whole cardiac cycle by adjusting the TD shift to the heartbeat. For each cardiac phase, three orthogonal diffusion directions at b-value 350 s/mm 2 were acquired. The other sequence parameters were kept the same as the ones used for the cDTI acquisition described above. The total scan time for the TD scout was~3 minutes for the corresponding 30 cardiac phases and 90 images.
In this work, three cardiac phases were studied: 1) early systole after the R-wave (TD~150-200 ms); 2) late systole as close as possible to maximum contraction (TD~250-350 ms); and 3) diastasis (TD~700-900 ms). Within these ranges, the TD for each cardiac phase was determined by visual assessment of the Trace signal (average across diffusion directions) of the TD scout acquisition. The TDs were identified by selecting the TD scout images with the least motion induced artifacts. Imaging all three cardiac phases in every volunteer was not always possible due to subject specific cardiac motion. However, at least two cardiac phases were successfully acquired in each volunteer.
For the retrospective quantitative analysis of the TD scout, the mean signal over the central zone of the DWI image was reported (see S1 Fig). Subsequently, the reported mean signals were normalized with respect to the maximum mean signal across all cardiac phases per volunteer.
For the cDTI acquisition, image post-processing included rigid registration, repetition averaging, and an in-plane interpolation using Fourier transform zero-filling. The twelve resulting diffusion encoded images and the non-diffusion weighted image were then fitted to the tensor model using a non-linear least-squares method [20]. The mean diffusivity (MD), the fractional anisotropy (FA), the eigenvectors (E 1 , E 2 , E 3 ), and the eigenvalues (λ 1 , λ 2 , λ 3 ) were extracted from the tensor model. The first eigenvector E 1 of the cDTI model represents the myofiber orientation, while the second eigenvector E 2 is usually attributed to the sheetlet direction. λ 1 , λ 2 , and λ 3 represent, the diffusivity along primary, secondary, and tertiary diffusion directions respectively.
Endocardial and epicardial LV regions of interest (ROIs) were manually drawn to exclude the blood pool, the papillary muscles, and the right ventricle (RV) (Fig 2A). The heart was divided into six segments defined by the American Heart Association (AHA) guidelines [21]. Voxel-wise helix angle (HA), E 2 angle (E2A), and transverse angle (TA) (Fig 2B) were computed after projecting E 1 onto the epicardial tangent (HA) and the horizontal plane normal to the epicardium (TA), and after projecting E 2 onto the plane normal to E 1 projection (E2A) [2,3,22] (Fig 2D). Lastly, a filter was applied to unwrap the HA (S2 Fig).

Quantitative analysis
Transmural HA distributions in the mid-ventricular slice for each subject were reported as a function of the normalized transmural wall depth measured from the epicardial and endocardial contours [22]. The HA distribution was described using nine points equally spaced along the transmural wall. For each point, the median HA was reported. For simplicity, only the epicardial, mid-myocardial, and endocardial medians are reported in the text and tables.
The HA distributions were also characterized by the HA range (HAR), which describes the angle difference between the endocardium and epicardium. HAR was calculated for each volunteer by subtracting the median HA at the epicardium from the median HA at the endocardium. Between two cardiac phases, the myofiber mobility was then calculated as the relative Finally, all quantitative values (ΔHAR, HAR, E2A, TA, MD, FA, λ 1 , λ 2, and λ 3 ) are reported as median [Quartile 1, Quartile 3] across all volunteers for the three cardiac phases. All statistical analyses were carried out in SPSS (Version 26, IBM Corp.) using a non-parametric Kruskal-Wallis one-way ANOVA. For statistically different distributions, corresponding pair-wise comparisons were realized between cardiac phases using a Bonferroni post hoc correction. p<0.05 was considered statistically significant.

Reproducibility of helix angle and E2A
An intra-and inter-observer study was performed to estimate the reproducibility of HAR and E2A. Two manual segmentations (intra-observer) were performed by two observers (interobserver) for a total of four observations per subject. Median HAR, E2A as well as intra-class correlation coefficient (ICC) were reported across subjects for the four observations. The first segmentation from the first observer was used for all non-reproducibility analyses reported in the results section.

TD scout and cDTI acquisitions
A total of nine healthy volunteers were scanned (4 males and 5 females, 28 ± 4 years old, 60 ± 11 beats per minute). First, the subject specific TDs were determined by using the TD scout sequence. Because they were during free-breathing (without a navigator) with less recovery time (TR~two heart beats instead of four), the TD scout images had lower signal and more non-motion related artifacts than the subsequent cDTI images. However, it was possible to visually identify motion related artifacts from the TD scout images. An example of the TD scout acquisition and the corresponding average signal per TD are displayed in Fig 3. As shown in Fig 3A and 3B, cardiac motion signal dropout was subject dependent and thus all three cardiac phases were not available for every volunteer. After visual assessment of the TD scout image quality, it was determined prospectively that cDTI acquisition was possible at late systole, early systole, and diastasis in 100%, 78% (7 out of 9 volunteers), and 67% (6 out of 9 volunteers) of the cases, respectively. The corresponding TDs across all volunteers for early systole, late systole and diastasis were 150 [148, 153] ms, 318 [293, 329] ms and 814 [800, 855] ms respectively. cDTI was more reliably acquired at late systole, while early systole and diastasis presented more motion artifacts ( Fig 3C). Four volunteers had high quality data during all three cardiac phases, three at early and late systole, and finally two at late systole and diastasis. The mean acquisition time for the high-resolution cDTI protocol per cardiac phase was 7 ± 2 min. Representative cardiomyocyte orientations E 1 , HA, E2A, and TA maps for one volunteer at all three cardiac phases are presented in Fig 4. Due to longitudinal cardiac motion, the three imaged cardiac phases may have been acquired at slightly different positions along the LV long axis. Using the DENSE acquisition from a subset of the volunteers (N = 5), the mid-ventricular longitudinal motion was estimated to be 4 mm from early to late systole, 4 mm from diastasis to early systole, and -8 mm between late systole and diastasis (positive values correspond to motion from base towards apex).

Helix angle and helix angle range
The HA distributions as a function of wall depth were extracted from the HA maps. Fig 5 shows an example HA distribution for the six AHA segments and for the entire mid-ventricular myocardium. A smaller IQR was observed in each AHA segment compared to the global HA distribution, demonstrating that the transmural HA distributions are segment specific. Additionally, the shape of the HA distribution for each AHA segment is largely conserved across each imaged cardiac phase as, for example, in the anterolateral segment (Fig 5). The epicardial myofibers in the anterior and inferior segments display a wide distribution due, in part, to the right ventricle (RV) insertions (Figs 5 and 6). Table 1 reports HA and HAR measured across volunteers in the six AHA segments. The intra-segment HAR was maximum at late systole compared to early systole and diastasis for the anterior, inferior, and anterolateral segments. In contrast, the intra-segment HAR was maximum in diastasis compared to early and late systole for the inferolateral and inferoseptal segments. ΔHAR between early and late systole was larger in magnitude (21.8 [17.1, 25.4]˚) than between late systole and diastasis (-4.2 [-1.3, 3.7]˚) or between diastasis and early systole (-11.8 [-29.0, 3.5]˚) as seen in Table 2. Overall, statistical differences were found in HAR across cardiac phases (p = 0.004). Pair-wise comparisons were significant between early and late systole (p = 0.005) and between diastasis and early systole (p = 0.004), but no statistical difference was detected between late systole and diastasis (p = 1). Overall, TA distributions were statistically different among cardiac phases (p = 0.033), but no pair-wise differences were observed after post-hoc correction between early and late systole (p = 0.094), late systole and diastasis (p = 1), and diastasis and early systole (p = 0.054).

Transverse angle and E2 angle
In contrast, an increase was observed in E2A distribution from early systole to late systole (    For E2A values, the ICC was 99% overall, across intra-and inter-observer analyses. Median values for E2A and HAR for the three cardiac phases per segmentation are given in S1 Table. Intra-and inter-observer variations did not affect the overall median nor the trends of E2A and HAR across cardiac phases.

Discussion
A second order motion compensated diffusion encoding gradient waveform design was used together with a SE approach to acquire cDTI images at a high spatial resolution of 1.6x1.6x8 mm 3 . To the best of our knowledge, this represents the finest spatial resolution acquired for in vivo human cDTI. The TD scout approach used in this work allowed a subject-specific calibration of the trigger delay for cDTI. Using this prospective calibration, two to three cardiac phases were identified for in vivo imaging of cardiomyocyte mobility. This approach characterized both global and regional cardiomyocyte orientations and demonstrated significant increases of HAR and E2A during contraction.

Multi-phase high resolution SE cDTI
The in vivo application of cDTI has been limited due to the sensitivity to cardiac bulk tissue motion, leading to the development of the STEAM sequence [23]. While effective at compensating for cardiac motion and providing a very good diffusion sensitivity, STEAM inherently leads to a significantly extended acquisition time, multiple breath-holds, and limited SNR-efficiency [9], all of which limit the spatial resolution [12]. The SE approach provides inherently higher SNR efficiency and can be used in conjunction with navigator-based free-breathing strategies that enable using a longer TR, thereby leading to more complete longitudinal magnetization recovery between acquisitions [24] for even more SNR. However, diffusion SE acquisitions require a motion compensated (M 1 = M 2 = 0) gradient waveform design, which incurs a TE penalty for a fixed b-value. This TE penalty can be partially reduced by optimizing the sequence parameters and the gradient waveform. In this study, the sequence parameters   1 [-9.6, 17.3]  9.7 [-13.0, 9.0]  -12.1 [-19.2, -4

PLOS ONE
were manually optimized to minimize the TE contribution from the EPI echo train length, which reduces the overall TE and EPI distortions. The diffusion encoding gradient waveform used a symmetric design to minimize M 2 based on the sequence parameters [13]. After this inline gradient waveform design, a residual dead-time (δ d , Fig 1C) of only 4ms remained. This dead-time could be further reduced by using a convex optimization approach for asymmetric gradient design [14]. Under these circumstances the benefits of using an asymmetric gradient waveform design would not have afforded significant TE savings. Gradient waveforms designed with a nulled second order moment are effective and efficient approaches to minimize artifacts induced by intravoxel constant velocity or constant acceleration. However, cardiac deformation is complex and differs from subject to subject. As a result, these linear motion assumptions do not apply to all cardiac phases nor all pixels, which can generate substantial diffusion signal loss in an unpredictable manner. The TD scout calibration was designed to efficiently identify cardiac phases with little signal dropout for each subject. This yields a fast tool for prospectively determining which cardiac phases are best suited for diffusion weighted imaging. As a prospective tool, the TD scout reduces exam time as it ensures that the subsequent cDTI acquisition is performed at the best cardiac phase, thereby avoiding long manual adjustments and scan iteration. Similar approaches to the TD scout have been used by Stoeck et al. [15] and Rapacchi et al. [25] to retrospectively identify acceptable trigger times for cardiac diffusion imaging.
The success rates obtained in this work (78%, 100% and 67% in early, late systole, and diastasis respectively) are close to the ones found by Scott et al [16] using a similar spin echo M1M2 motion compensation as the one used in the current study. The success rate below 100% in early systole and diastasis shows the limitations of this motion compensation approach. As shown in Fig 3, signal dropout due to cardiac motion was observed, particularly in diastole while most of the systolic phases were available for imaging. The cardiac phase imaged in diastole corresponded to diastasis, where the LV motion pauses during filing before atrial systole begins. However, diastasis usually shortens as the heart rate increases and the beat-to-beat variation in cardiac motion may not respect the linear velocity/acceleration assumption behind the M1M2 waveform. Previous works have shown that non-motion compensated waveforms combined with a signal recovery algorithm [25,26] were able to reliably scan in diastole because of a significantly shorter TE leading to a very small temporal footprint. Typical in-plane spatial resolutions in clinical cDTI studies range from 2x2 mm 2 to 3x3 mm 2 , usually limited by low SNR. Because HA varies transmurally, higher spatial resolution could play an important role in improving the assessment of regional cardiomyocyte orientation and mobility. Previously McClymont et al. [27] have demonstrated that typical spatial resolutions could underestimate the HAR by up to 18˚, in part due to partial volume effects. In this study, high spatial resolution enabled the detailed analysis of in vivo HA per AHA segment, revealing distinct microstructural organization that varies per-segment. However high resolution remains challenging, particularly in the inferolateral segment which is prone to strong geometrical distortion due to the lung/heart/liver interface [28].

Microstructural mobility
In this work, both global and regional HA mobility and a significant increase in HAR during contraction were observed. The HAR in early systole 76.9 [72.6, 80.5]˚, late systole 96.6 [85.9, 100.3]˚, and in diastasis 91.7 [85.9, 100.8]˚agree with early reports of ex vivo cDTI animal data. Omann et al. [29], using submillimeter ex vivo DTI in swine hearts, reported an LV HAR of 80˚and 71˚in contracted and relaxed states. These HAR trends are repeatedly observed in other ex vivo swine studies, with reported HAR of 109˚and 89˚ [8], and 107˚and 87˚ [30] in the contracted and relaxed states, respectively. In rat hearts, Chen et al. [31] reported a HAR of 111˚, 134˚, and 105˚in early systole, peak systole, and diastole, respectively. In vivo cDTI studies of the healthy human LV have also shown a steepening of HA during systole. Stoeck et al. [22] reported a HAR of 56˚and 72˚in diastole and systole, while von Deuster et al. [32] reported a HAR of 54˚and 77˚.
A detailed characterization of healthy HA mobility is valuable in providing a baseline to study the microstructural mechanisms of LV dysfunction. Indeed, several studies have implicated deleterious HA remodeling and abnormal mobility in heart disease. For example, an increase in the transmural slope of cardiomyocyte elevation was recently described in infarcted tissue [33] with respect to remote myocardium. von Deuster et al. [32] observed a decrease in HA mobility from diastole to systole in patients with dilated cardiomyopathy. Recently, Gotschy et al. [11] demonstrated the link between global longitudinal strain and HA mobility in patients with cardiac amyloidosis. Clinical studies on hypertrophic cardiomyopathy [3] appear to show a steepening of HA during contraction, although the reporting primarily focused on E2A changes.
We also observed an increase of E2A in late systole compared to early systole and diastasis. As shown in Fig 2, E2A represents a combination of projections involving both E 1 and E 2 , and is usually attributed to the cardiomyocyte sheetlet direction [3]. In this work, E2A increased from 27.7 [20.8, 35.1] 20.7 [16.6, 26.4]˚in diastasis. Using a STEAM-based approach and at slightly different cardiac phases, studies in healthy controls demonstrate E2A changing from 65˚to 18˚ [8] and from 56˚to 24˚ [3] from peak systole to diastole, respectively. In the same studies and in patients with hypertrophic cardiomyopathy, E2A was reported to decrease from 74˚to 48˚ [8] and from 64˚to 47˚ [3] from peak systole to diastole, respectively. E2A changes were also observed in patients with dilated cardiomyopathies [8] and amyloidosis [11], showing the importance of measuring E2A in addition to HA as a diagnostic biomarker and to better understand how changes in E2A relate to ventricular contraction.
Although HAR and E2A mobility trends are consistent across cDTI studies, there is a broad range in the reported magnitude of HAR and E2A changes during contraction and filling. As shown by Scott et al [16], these variations could be attributed in part to difference in sensitivity (diffusion time) between SE and STEAM imaging. In addition, a collection of artifacts (motion, off-resonance, chemical shift of fat, etc.) are also present in cDTI, with a different effect in SE and STEAM imaging [9,16,22,34]. Finally, the angular metrics used to represent the microstructure (namely HA, TA, E2A, HAR) are based on circumferential and radial directions constructed from manually defined segmentations. This inter-observer dependence is particularly significant in evaluating the HAR metric based on the HA in proximity of the endocardial and epicardial surfaces. This can be seen in the intra-and inter-observer study presented in S4 Fig, where we report more variation across subjects for HAR than E2A. However as shown in S1 Table, median trends during the cardiac cycle remain consistent for both intra and inter-observer analyses.
Cardiomyocyte mobility has been reported in a wide range of cardiac phases across cDTI literature. As shown in the TD scout measures SE approaches can measure most of contraction up to late systole (~150 to 300ms), and during the diastasis phase in diastole. In this work, the complete cardiomyocyte mobility cycle was not observed, as the cardiomyocytes in diastasis do not completely return to their initial configurations in early systole. The return in cardiomyocyte orientations to complete the cardiac cycle must occur outside the window of the cardiac phases captured in this work. Future studies are needed to determine precisely when HA and E2A return to their configuration at the beginning of systole.

Diffusivity values at early systole, late systole, and diastasis
The MD and FA values reported in this work are in the range of previously reported values in healthy volunteers using spin-echo sequences.  [22], a higher MD was found in systole 1.2 ± 0.2 x 10 −3 mm 2 /s compared to diastole 0.9 ± 0.1 x 10 −3 mm 2 /s using a STEAM approach. Scott et al. [16].found the opposite trend with an MD of 1.02 [0.14] x 10 −3 mm 2 /s in systole, 1.07 [0.12] x 10 −3 mm 2 /s at the systolic sweet spot and 1.13 [0.08] x 10 −3 mm 2 /s in diastole using a STEAM approach. In the same study using a SE sequence, MD was 1.46 [0.43] x 10 −3 mm 2 /s in systole, 1.70 [0.18] x 10 −3 mm 2 /s at the systolic sweet spot and 1.78 [0.34] x 10 −3 mm 2 /s in diastole. In our previous work [35] using SE, a higher MD was found in diastole 1.91 ± 0.34 x 10 −3 mm 2 /s compared to systole 1.58 ± 0.09 x 10 −3 mm 2 /s. Values were reported as mean ± standard deviation, Median [Interquartile Range] or Median [Quartile 1, Quartile 3]. The source of this increase could be explained by several effects. The higher MD in early systole could be partially due to micro-perfusion, which is known to increase during filling and whose signal could not be fully spoiled by the motion compensated diffusion encoding gradients [36]. Despite our use of a motion compensated waveform, unwanted cardiac bulk motion present in early systole could artificially increase MD [15]. However, in this work, the diffusivity differences are unlikely due to a change in cellular/extra-cellular compartments between cardiac phases. The b-value (350 s/ mm 2 ) and the diffusion encoding time used here (Δ = 40ms) were too small to sensitively capture these compartmental changes [16], which is also confirmed by fact that the ratio λ 2/1 remains constant across cardiac phases.

Limitations
A first limitation of this study is the single-slice heart coverage acquired with a long TR to maximize SNR. This approach limits characterizing cardiomyocyte mobility throughout the heart and limits overall SNR efficiency. Several approaches have been proposed to improve the spatial coverage without increasing scan time, such as simultaneous multi-slice acquisition [37] and the slice-following navigator technique [24]. Both temporal and spatial coverage could be improved using advanced dual cardiac phase encoding [22]. All these approaches would improve SNR efficiency and could offer improved clinical applicability of the presented protocol.
Second, all three cardiac phases of the cDTI datasets were acquired at the same spatial location for each volunteer. However, the myocardium passes through the imaging plane as the heart contracts and fills and thus the tissue imaged in the three cardiac phases may differ slightly. The mid-ventricular longitudinal motion from base to apex was approximately 4mm from early to late systole, 4mm from diastasis to early systole and -8mm between late systole and diastasis (estimated using DENSE imaging on a subset of the study subjects). Considering the slice thickness (8mm) and the relative longitudinal homogeneity of the cardiomyocyte orientation [2], it is possible that through-plane tissue motion mildly influences the microstructural mobility reported in this cDTI study. Few multi-cardiac phase in vivo cDTI studies have considered the longitudinal displacement of the heart by estimating this displacement prospectively [3].

Conclusion
We demonstrate that it is possible to acquire high-resolution cDTI datasets in~7 minutes per cardiac phase during free-breathing. This is compatible with clinical time constraints and suitable for patients with respiratory problems. In healthy volunteers, in vivo regional HA mobility was reported at three distinct cardiac phases and suggests that on average cardiomyocytes reorient more longitudinally during contraction, while sheetlets tilt in the direction of LV wall thickening. In the future, these observations may provide clinical insight into LV function in healthy and diseased hearts. Given the intricate structural and microstructural phenotypes associated with heart disease [4,5,38], high resolution in vivo cDTI has the potential to offer a better understanding of the complex link between microstructural remodeling and the underlying cardiac dysfunction at a patient-specific level.  -observer). The intraclass correlation coefficient (ICC) for HAR was 76% overall and 84% inter-observer. The intra-observer ICC was 72% for Observer-1 and 92% for Observer-2. The ICC for E2A was 99% overall, inter and intra-observer. (TIF) S1 Table. Median Helix Angle Range (HAR) and E2 Angle (E2A) computed by two observers across volunteers. Each observer segmented the data twice. Intra-and inter-observer variability does not affect the overall changes in HAR and E2A across the analyzed cardiac phases. (TIF)