Imaging the dynamics of cardiac fiber orientation in vivo using 3D Ultrasound Backscatter Tensor Imaging

The assessment of myocardial fiber disarray is of major interest for the study of the progression of myocardial disease. However, time-resolved imaging of the myocardial structure remains unavailable in clinical practice. In this study, we introduce 3D Backscatter Tensor Imaging (3D-BTI), an entirely novel ultrasound-based imaging technique that can map the myocardial fibers orientation and its dynamics with a temporal resolution of 10 ms during a single cardiac cycle, non-invasively and in vivo in entire volumes. 3D-BTI is based on ultrafast volumetric ultrasound acquisitions, which are used to quantify the spatial coherence of backscattered echoes at each point of the volume. The capability of 3D-BTI to map the fibers orientation was evaluated in vitro in 5 myocardial samples. The helicoidal transmural variation of fiber angles was in good agreement with the one obtained by histological analysis. 3D-BTI was then performed to map the fiber orientation dynamics in vivo in the beating heart of an open-chest sheep at a volume rate of 90 volumes/s. Finally, the clinical feasibility of 3D-BTI was shown on a healthy volunteer. These initial results indicate that 3D-BTI could become a fully non-invasive technique to assess myocardial disarray at the bedside of patients.

The myocardial fibers architecture plays a major role in the cardiac function. Fibers orientations vary continuously and smoothly through myocardial walls 1 and this complex organization is closely linked to the mechanical and electrical myocardial function [2][3][4] . For instance, the helicoidal fiber distribution in the left ventricular walls contributes to the torsional motion of the heart during the ejection phase, which is believed to improve the heart's pumping efficiency 5,6 . The fibers of the heart are also related to its electromechanical properties [7][8][9][10] ; indeed, the electrical activation propagates preferentially in the direction of myocardial fibers and results in their synchronous contraction. Myocardial fiber disarray is thought to appear in the early stage of many pathologies such as in cardiomyopathies or in fibrosis 11 . The mapping of the fiber architecture, on the one hand, could increase our knowledge of the cardiac function 12 , and on the other hand, could potentially enable early diagnosis of cardiomyopathies. Yet, no method to map the fiber architecture is currently used on a regular basis for clinical purposes.
Magnetic Resonance Diffusion Tensor Imaging (MR-DTI) 13,14 is widely used for mapping the structural connectivity of the human brain. However, its application to the heart, while feasible in animals 15 and humans 16 , remains challenging in vivo due to limited frame rates (and hence long acquisition times) and limited robustness of the MR-DTI signals to tissue motion 11,15,16 . Optical methods such as optical coherence tomography 17 and two-photon microtomy 18 can also map the fiber directions at the microscopic level but remain limited to superficial ex vivo tissue and to small regions of interest. In ultrasound imaging, several methods have been proposed to quantify the anisotropy of various physical parameters linked to the fiber orientation, including the ultrasonic attenuation 19,20 , the integrated backscattered intensity [19][20][21][22] and the myocardial stiffness using Elastic Tensor Imaging (ETI) 23 . However, none of these methods has been yet implemented for the time-resolved, 3D mapping of the myocardial fibers orientation.
In this paper, an entirely novel ultrasound-based imaging technique called Backscatter Tensor Imaging (3D-BTI) is described and evaluated in vivo with respect to its capability of mapping the myocardial fibers orientation and its time-resolved dynamics during an entire cardiac cycle. 3D-BTI is based on ultrafast volumetric ultrasound acquisitions, which are used to quantify the spatial coherence of backscattered echoes at each point of the imaged volume at a rate of 90 volumes/s. Spatial coherence is a physical property of backscattered waves that depends on the distribution of scatterers at the subwavelength scale. Ultrasonic spatial coherence has thus the capability to provide statistical information on the tissue microstructure at a scale far below the resolution of ultrasound images. In 3D-BTI, the anisotropy of spatial coherence is analyzed on a 2D matrix probe in order to derive its principal directions and derive the fibers orientation. The link between the anisotropy of spatial coherence and the fibers orientation was previously demonstrated in composite materials and in fibrous soft tissues such as the myocardium and the skeletal muscle 24,25 .
Another key aspect of 3D-BTI is the use of volumetric ultrafast plane wave imaging. We recently developed a 3D fully programmable ultrasound system for volumetric ultrafast plane wave imaging achieving volume rates of up to 5000 volumes/s 26 . By coherent compounding of the backscattered echoes associated with tilted plane wave emissions, the spatial coherence and thus the fibers orientation can be obtained simultaneously in each voxel of the volume at high volume rate.
3D-BTI is hereby found to be capable of mapping the transmural fibers orientation in a beating heart. For the purpose of this validation study, 3D-BTI was compared to histological analysis on ex vivo explanted myocardial tissues. An excellent agreement on the transmural fiber orientation was found between 3D-BTI and histology. The in vivo feasibility of 3D-BTI on a beating heart was then demonstrated in an open chest ovine model. Finally, transthoracic imaging was performed on a healthy volunteer to show the clinical feasibility of 3D-BTI. To our knowledge, it is the first time that the dynamics of the myocardial fibers are observed non-invasively over an entire cardiac cycle of a single heartbeat in humans.

Results
The general principle of 3D-BTI is illustrated in Fig. 1. First, tilted plane waves (Fig. 1A) are emitted from a 2D matrix array probe connected to a customized, programmable, ultrasound system 26 . For each emitted plane wave, the backscattered echoes received by each element of the matrix probe are recorded and further processed using 3D plane-wave coherent compounding to synthetically generate a voxel-specific focal region (Fig. 1B). The spatial coherence associated to each focal region is then computed at each point of the 3D volume (Fig. 1C) and used to determine the fiber orientation by applying an elliptic fit (Fig. 1D). Finally, a vector representation is used for the visualization of fibers in 3D space (i.e. Fig. 1E).

Ex vivo experiments.
To establish the accuracy and precision of 3D-BTI, myocardial samples explanted from the anterior wall of the left ventricle of a porcine heart were imaged using 3D-BTI and further analyzed by histology. Figure 2A shows the 3D grayscale ultrasound image of one myocardial sample and three examples of spatial coherence functions obtained at different depths ( Fig. 2B,C,D). At each location, the spatial coherence presents a strong anisotropy with a principal direction that varies with depth. The vector representation (Fig. 2F) of the entire volume shows that fiber orientation varies across the wall. Fiber orientation was found to vary gradually through the ventricular wall, which is consistent with the literature 1 . The transmural orientation was averaged over the volume for the 5 myocardial samples and was found to vary continuously through the wall with an average difference of 98.6° ± 8.9° between endocardium and epicardium. These results were compared against the fiber orientation obtained from histological analysis. Figure 3A and B shows a representative sample studied with both histology and 3D-BTI. Figure 3C shows a Bland-Altman analysis applied to all five samples of this study. A bias of 0.33 degrees and 95-% limits of agreement equal to −10.1° and 10.7° were found.

In vivo open chest experiments. The in vivo feasibility of 3D-BTI was demonstrated in the beating heart
of an open chest sheep. 3D-BTI acquisitions were performed during a complete single cardiac cycle with a volume rate of 90 volumes/s. Figure 4 shows examples of fiber reconstruction at specific times of the cardiac cycle (A Late Diastole, B Early Systole, C Late Systole, D Early Diastole). It demonstrates the feasibility of following the fiber orientation during an entire cardiac cycle at high frame rate (the complete cine-loop is shown in Supplemental Movie 1). The transmural fiber distribution at one location as a function of time is shown in Fig. 5. Fiber angles were found to vary as a function of depth with an absolute difference of 96° from the epicardium to the endocardium.
In vivo transthoracic BTI of the human heart. Finally, the feasibility of 3D-BTI on the human heart was shown on a healthy volunteer in a realistic clinical transthoracic imaging setting. Transthoracic conventional gray-scale imaging was used to position the probe on the parasternal view. 3D-BTI was then performed with a total imaging depth of 60 mm. Figure 6 shows the helicoidal fiber distribution imaged in the left ventricle anterior wall at mid-level both in systole and in diastole. Fiber angles were found to vary transmurally with an absolute difference of 104 ± 9° from the epicardium to the endocardium.

Discussion
In this study, we have introduced 3D-BTI, a novel ultrasound-based technique for the mapping of the fibers orientation in myocardial tissues. The objectives were to determine the accuracy and precision of 3D-BTI in mapping the fibers structure and to demonstrate its feasibility on beating hearts. The accuracy and precision of 3D-BTI were validated in 5 ex vivo porcine left ventricular myocardial samples. The fibers angles distribution of each sample was obtained within volumes of 3 cm 3 and the transmural angle variations of each sample were compared to histology. An excellent agreement was found against histology, therefore demonstrating the accuracy and precision of 3D-BTI.
The in vivo feasibility of 3D-BTI was then demonstrated in a beating sheep heart. The fibers distribution of the left ventricular free wall was successfully imaged during one entire cardiac cycle at a volume rate of 90 volumes/s. To our knowledge, this is the first time that the dynamics of the myocardial fibers were mapped during an entire cardiac cycle at high frame rate. The transmural fiber variation was found to be 96° in diastole and did not change significantly during the systolic phase despite a 56% increase in wall thickness. This small variation over the cardiac cycle confirms the results of several studies that have investigated the fibers orientation in diastole and in systole using histology 27 . Finally, the feasibility of BTI in a realistic clinical transthoracic imaging setup was shown on a healthy volunteer and allowed for the mapping of the myocardial fibers orientation of the human heart both in systole and in diastole.
3D-BTI could not only bring new insights in the knowledge of the cardiac mechanics but it could also become a major tool to non-invasively detect fiber disorders linked to cardiac diseases in early stages. Myocardial fiber disarray has been shown in post-infarct remodeling of the left ventricle as well as in hypertrophic cardiomyopathy. 3D-BTI may be used to quantify disarray in the fiber angles but also to map the fractional anisotropy in order to quantify collagen infiltration and fibrosis content. Fractional anisotropy provides a measure of the level of anisotropy and can be obtained using the coherence values along and across the fibers.
3D-BTI presents a number of advantages as it can be performed non-invasively at high frame rate, and in vivo in beating hearts. In contrast to shear wave imaging based techniques, 3D-BTI is based on low energy emissions and thus can be performed continuously to provide the complete description of the fiber dynamics at high frame rate. Moreover, while the present study focused on the use of plane waves for the generation of synthetic foci, diverging waves 28,29 , which can achieve even larger fields of view and image the entire heart in 3D 26 , could also be used to improve 3D-BTI. The use of diverging waves would, however, result in a non-uniform spatial resolution. This is a limitation that needs to be evaluated in more detail.
3D-BTI is robust to motion as a result of the use of ultrafast volumetric acquisitions, which enabled the mapping of the coherence functions in entire 3D volumes using a limited number of transmits. For example, while 49 plane wave emissions were sufficient to perform 3D-BTI in each voxel, performing the equivalent acquisition using standard focusing would have required 64*64*500 emissions, which would corresponds to a few minutes per volume instead of a few milliseconds in the case of plane wave imaging. One of the current limitations of 3D-BTI is its incapability to map the z component of the fiber direction. Despite this limitation, mapping of fiber directions remains possible in anisotropic soft tissues in which the fibers are parallel to the 2D array plane, which is the case for the skeletal muscle or the anterior wall of the myocardium in standard parasternal short and long axis views. However, to apply 3D-BTI to other tissues in which structures can be oriented along the z-axis in sagittal and coronal views, tilted synthetic foci would be required. This could be achieved using subapertures in reception 30 , which is the object of on-going work in our group. Moreover, while the ultrasound acquisitions required to perform 3D-BTI are very short (~10 ms), postprocessing operations required up to twenty minutes to achieve the vectorial representation of a full volume. Further optimization of the algorithms along with the rapid growth of computational power of video cards and processors could potentially lead to a real-time implementation of 3D-BTI in the future. Segmentation of the epicardium and endocardium remained challenging using the ultrafast Bmode images and required to be performed manually by a trained  BTI may be affected by low ultrasound signal to noise ratio at large depths. The performance of BTI need to be investigated in more details for low SNR configurations. Nevertheless, we can anticipate BTI to be robust to incoherent noise such as clutter and thermal noise due to the suppression of a large part of incoherent noises by the coherence function. However, we can also anticipate that the performance of BTI will drop considerably in very low SNR situations where incoherent noise becomes larger that the coherent signals.
Finally, while presenting similar advantages in terms of portability, safety, and real-time imaging capabilities, the 3D ultrafast ultrasound system that we used is a unique laboratory prototype that is not currently available in clinical practice. In summary, 3D-BTI may constitute a unique tool for non-invasively evaluating the myocardial fiber structure of a patient suffering from myocardial fibrosis or hypertrophic cardiomyopathy for diagnosis, treatment monitoring and follow-up.

Materials and Methods
Spatial coherence on a matrix array transducer. Ultrasonic spatial coherence was assessed experimentally by focusing an ultrasound wave in a medium and by receiving the backscattered echoes on a 2D matrix array transducer (Fig. 7A and B). The 2D spatial coherence function R(Δx, Δy) was obtained by computing the auto-correlation of the signals received by pairs of elements i and j distant by Δx and Δy, along the main two coordinate axis of the matrix array: where S k (x k , y k , t) is the signal received on the element k of the matrix transducer with coordinates x k , y k after applying a time-delay to compensate the propagation path. Δx and Δy are the distances between and element i and element j: (2) i j i j T 1 and T 2 represents the averaging temporal window (equivalent to depth). In random media, i.e. in absence of fibers, ultrasound scatterers are randomly distributed, (Fig. 7A), so that the coherence function is predicted by the so-called Van Cittert-Zernike theorem 31 : the coherence function is given by the spatial Fourier transform of the intensity distribution at the focal spot. For a square matrix array, the pressure distribution at the focal spot is a 2D sinc function so that the coherence function appears as the Fourier Transform of a squared 2D sinc function. The pyramidal shape coherence function of a N x × N y 2D matrix array transducer is displayed on Fig. 7C.
Data acquisition and signal processing. In order to focus the ultrasonic wave at specific locations in the medium, conventional ultrasound imaging relies on applying time delays to the transmitted wave (Fig. 7). A large number of transmitted waves focused at different locations of the medium is then required to obtain one single image of the medium. In volumetric imaging, the number of transmitted waves can reach several thousands of waves which results in very low volume rates. To overcome this issue, our approach consisted in using ultrafast plane wave imaging with a coherent compounding approach 26 . Coherent compounding can generate, at any location of the region of interest, synthetic focal zones using only a few tens of emissions. Therefore, with this approach, large volumes can be imaged at very high volume rate between 50 and 5000 volumes/s. The emissions consisted in the transmission of several tilted 2D plane waves (i.e. Fig. 1A). Each plane wave was defined by two angles. All pairs of angles ranging from −6° to 6° with a step angle of 2° were emitted, for a total of 49 tilted plane waves. Plane waves were emitted from a 2D matrix array probe with 32 (x-axis) by 35 (y-axis) elements (3 MHz, 0.3 mm pitch, 50% bandwidth at −3 dB, Vermon) connected to a customized, programmable, 1024 channel ultrasound system described in Provost et al. 26 . This system was designed and built in-house for 3D ultrafast imaging applications. This system allowed us to perform the first in vivo 3D ultrafast imaging acquisitions in various applications such as 3D imaging of the cardiac blood flows 26 , 3D elastography of soft tissues 32 and 3D mapping of blood vessels 33 . The 1024 independent channels could be used simultaneously in transmission, whereas receive channels were multiplexed to 1 of 2 transducer elements. Therefore, each emission was repeated twice, with the first half of the elements receiving during the first emission, and the second half of the elements receiving during the second emission.
The received signals were sampled at 12 MHz, for all the 1024 elements and the radio-frequency data were recorded. The 3D images were computed off-line using a dynamic receive focusing beamforming algorithm followed by coherent compounding (see Fig. 1B). The 3D volume had a lateral size of 9.6 × 9.6 mm and a 3 cm to 6 cm depth depending on the region of interest. The lateral sampling of the image was 0.3 mm laterally (32 × 32 lines) and the axial sampling was 0.05 mm.
2D coherence functions were computed at each point of the 3D volume (Fig. 1C) using equation (1). The temporal window [T 1 T 2 ], was set to 1.6 µs (i.e. five periods at 3 MHz). The coherence functions were averaged at each depth of the 3D volume using a lateral sliding window with a kernel of (0.15 × 0.15 mm). The algorithm was implemented on parallelized Nvidia GPU using CUDA language.
An algorithm based on the Radon transform was applied on the 2D coherence functions, to determine the fiber orientation. The first step of the algorithm consisted in computing the integral of the 2D coherence function along a line y = 0. The line was then rotated step by steps in the 2D coherence function plane, to compute the integrals at each angle in the range 0-179 degrees. The fiber orientation was determined by the angle that corresponded to the maximum value of the integral (i.e. Fig. 1D). In this study, steps 1 degree were chosen. Finally, a vector representation was used to visualize the fibers directions in 3D using the Amira software (Visualization Sciences Group, Burlington, MA) (i.e. Fig. 1E).

Experimental Setup. 1) In vitro experiments and histology.
In vitro experiments were performed on five porcine fresh left ventricle myocardial samples embedded in a 2% gelatin gel. A cross-marker whose orientation matched the local cardiac circumferential-longitudinal coordinates was labeled on the myocardial sample. The probe was positioned at 15 mm of the epicardium. Each myocardial region examined by BTI was dissected from the intact porcine heart into a rectangular block (10 × 10 × 25 mm). This block was fixed in formalin for 48 h. To keep the orientation of the tissue block in the wright direction, the basal aspect of the tissue block was labeled with tattoo ink. Then the tissue block was cut in 4 to 5 tissue slices which were embedded in paraffin. Each paraffin block was completely cut at 5 µm of thickness in serial sections. One histological section out of 50 was stained with H&E: the resulting gap between two stained adjacent sections was 250 µm. All the sections were scanned with a Hamamatsu Nano Zoomer (Hamamatsu City, Japan). Fiber angles in all the histological digitalized images studied were computed using the Hough transform 34 .

2) In vivo experiments.
Animal experiments were performed in an open chest sheep on the mid-anterior region. The experimental procedure was approved prior to use by the Institutional Animal Care and Use committee of Hôpital Européen Georges Pompidou (Paris Descartes) according to the European Commission guiding principles (2010/63/EU). 3D-BTI was thus performed within an entire cardiac cycle. The 3D-BTI acquisition was performed at a volume rate of 90 volumes/s during 1.2 seconds in order to cover more than one cardiac cycle. It enabled the reconstruction of 100 fiber volumes. During the experiment, the electrocardiogram was recorded.
Human experiments were performed by a trained cardiologist. The volunteer was positioned on the left lateral decubitus position. The probe was placed in the parasternal view and 2D real-time imaging was performed to position the probe and image the antero-septal wall. The 3D BTI acquisition was then launch to acquire 1 s of data at a rate 90 volumes/s. The study was carried out in accordance with the Declaration of Helsinki and was approved by the French ethical committee (CPP Paris Ile de France VI, N° 15-15). Informed consent was signed by the volunteer.