1. Introduction

The heart is the first organ to function in vertebrates during development. As a pump, its beating dynamics and the resulting mechanical outputs, such as the contractile and flow-induced forces, are essential for proper development of the cardiovascular system and the entire embryo [1–3]. Studying the cardiodynamics and hemodynamics in the embryos of animal models is therefore important to understand cardiogenesis and congenital heart defects, where live imaging plays a critical role in revolutionary discoveries [4–9]. Particularly, due to the dynamic nature of heart, fast imaging able to reconstruct the beating embryonic heart in 3D has been vital to reveal its structure, morphology and function, accompanying various molecular genetic tools to enable innovative studies in the early heart development.

Multiple high-resolution, high-speed optical modalities are utilized for dynamic imaging of the embryonic heart [10–17], with advantages and applications for specific model organisms. Among these, optical coherence tomography (OCT) presents a unique imaging scale and contrast that complement the widely used approaches, such as confocal and light-sheet microscopy. Using near-infrared light and low-coherence gating, OCT probes the optical scattering of tissue for label-free imaging and provides a micro-scale resolution with a millimeter-level depth [18]. This scale makes it possible for 4D (3D + time) visualization of the detailed structure of the early heart in live embryos [19–22], especially in the mouse model where imaging embryonic cardiodynamics within the intact yolk sac is difficult with other methods. With Doppler OCT, 4D intracardiac blood flows can also be imaged and quantified together with the structure information [23–26]. Notably, this has enabled spatiotemporal assessment of the flow-induced shear stress [27,28], studying the role of blood flow patterns in cardiac malformations [29,30], and investigating the pumping function of the valveless embryonic heart [31,32]. Such dual-contrast, high-resolution cardiodynamic and hemodynamic imaging with OCT is powerful for dynamic phenotyping and biomechanical analysis of the early cardiogenesis.

The imaging speed of OCT is important for capturing the beating embryonic heart. During early development, the heart rate generally falls in the range of 1-3 Hz. However, a sampling frequency well beyond the Nyquist frequency of this range is required to reveal the features of the heartbeat. For example, in the study with confocal imaging of zebrafish embryos [4], a 151 Hz sampling rate (more than 50 times the heart rate) allowed the researchers to identify the bi-directional, high-speed contraction wave that is foundational in discovering the new suction mechanism of valveless pumping. This type of wave traveling within a ∼200 µm distance would not be detectable under a much slower imaging rate. Using OCT, repeated 3D acquisition of the heartbeat represents a direct approach for 4D imaging [22,33,34]. While the continuous improvement of the OCT imaging speed has led to a direct volumetric embryonic heart imaging at ∼43 Hz [34], due to the interdependence between the rates of spatial and temporal sampling, the total amount of voxels remains limited even with a MHz-level A-scan rate [34]. This wide sampling interval makes further, in-depth analyses challenging.

In contrast, post-acquisition synchronization of the heartbeat phase constitutes another approach for 4D OCT imaging of the embryonic heart in which the spatial sampling rate can be largely decoupled from the volume rate [35–38], permitting high resolution in both space and time. The approach relies on an image-based alignment of the time-axis for repeated OCT B-scans acquired at different locations of the beating heart. The scan protocols include the generic B-M scans [35,36] along with the rotational scan of time-lapse B-scans [37,38], where the gating signal/image comes from the B-scan or M-mode cross-correlation (similarity) [35,36] in the case of the former, and from the common central A-scan [37,38] in the case of the later. As a result, the 4D volume rate equals the B-scan imaging rate, which can easily reach 100 Hz or more, and the spatial sampling is minimally affected by the imaging speed. This high sampling rate makes it simple to incorporate Doppler OCT that requires specific imaging parameters [26], such as a dense spatial sampling of A-scans. In comparison with ultrafast, direct volumetric 4D imaging, the post-acquisition synchronization approach provides a superior imaging quality, which proved sufficient to study the heart [29,31,32,39,40]. Moreover, it lowers the specifications for the OCT system without the need for high-end, ultrahigh-speed systems.

Recently, consumer-level OCT systems have become increasingly accessible to a wide range of biomedical labs and researchers, thanks to the efforts in commercialization and low-cost designs [41,42]. However, to enable a broader use of OCT as a 4D imaging tool to study the embryonic heart, especially among biologists, a post-acquisition synchronization method that is more accessible, efficient, and accommodating of the simplest OCT scan scheme is required.

Here, we report an open-source, post-acquisition synchronization method with OCT for 4D cardiodynamic and hemodynamic imaging of the embryonic mouse heart. Relying on the difference between images sampled adjacently in time as a measure of heart wall movement, this method provides good sensitivity to the cardiac activity when aligning heartbeat phases, even at early stages when the heart wall is thin, occupying only a small number of pixels. The method works with a single linearly scanned 3D data acquisition, in contrast with repeated B-scans at each scan location, the former of which is readily available in most commercial systems. We present the reconstruction quality of the new method as a comparison with the existing method [43] that has previously been demonstrated for OCT imaging in the mouse model [21,26,44]. We show the robustness of the method through 4D imaging of different mouse embryos with heart rates ranging from 1.24 Hz to 2.13 Hz. With the alignment process conducted on a 1D profile, the method has high efficiency, featuring sub-second alignment time while utilizing ∼100% of the original image files. This speed allows us to reconstruct repeated, dual-contrast 4D imaging of the beating mouse embryonic heart over a two-hour development period. We designed and implemented this method to maximize its accessibility, which will facilitate more research employing the advantages of OCT in the study of early cardiogenesis.

2. Materials and methods

3.1 Mouse embryos

We set up timed mating of CD-1 mice (Charles River Laboratories) and checked the vaginal plug in the morning. The day of observing the plug was denoted as embryonic day 0.5 (E0.5). We performed dissection and static culture of live embryos within the intact yolk sac on ∼E8.5-E9.0 following the established protocols [45,46]. Embryos were dissected with a small portion of the decidua left attached to the yolk sac, allowing us to position the embryo and keep the entire sample stable. To demonstrate the 4D reconstruction method for different heart rates, we utilized four embryos at developmental stages spanning ∼E8.5-E9.0. All animal manipulations were approved by the Institutional Animal Care and Use Committee at Baylor College of Medicine, and the experiments followed the approved procedures and guidelines.

3.2 OCT system and imaging

We performed embryonic imaging using a customer-built, spectral-domain OCT system with a central wavelength of ∼810 nm [47]. The system had an axial resolution of ∼5 µm in tissue and a traverse resolution of ∼4 µm. We placed the sample arm of the system inside an incubator with 37 °C and 5% CO₂ for imaging the embryo culture. As shown in Fig. 1(A), we positioned the embryo with the heart facing up for upright imaging through the culture medium. The imaging scale of OCT is ideal for mouse embryos at the stages both before and during cardiac looping. The imaging depth covered the entire embryonic heart in the intact yolk sac, including the whole embryo at the early stage (Fig. 1(B)). Without needing to open the yolk sac, the experimental procedure was free of extensive embryo manipulations, and prolonged embryonic development can be achieved in static culture for up to 24 hours [48].

 

Fig. 1. (A) An illustration of the OCT sample arm inside an incubator for imaging the beating heart in the live mouse embryo within the intact yolk sac. (B) An OCT depth-resolved image of an E9.0 live mouse embryo within the intact yolk sac. The scale bar is 200 µm.

Download Full Size | PDF

The new 4D synchronization method was designed for a dataset acquired over the beating heart as a single 3D image. The acquisition requirement is for B-scans within the 3D image to be densely sampled with a small spatial step, and for the total acquisition time to be long enough so that multiple heartbeat cycles are captured while the B-scan plane crosses every step within the transverse resolution of OCT. In this study, we performed OCT imaging of the beating heart at 14.7 µs per A-scan with a B-scan rate of either 100 Hz or ∼83.3 Hz. We maximized the spatial sampling density within B-scans for detecting Doppler OCT signals from the blood flow. In the data acquisition for 4D reconstruction, we acquired a single 3D image with a linear scan of the B-scan location that has dense spatial sampling in the direction across the B-scan. As shown in Fig. 2(A), the B-scans were acquired in the third space axis, Z, but also in time, t, due to the slow progression of the Z galvanometer. In this paper, we used only “Z” to represent this spatiotemporal axis. The Z step size (between B-scans) was selected such that the distance of the imaging beam movement over one heartbeat cycle was less than or equal to the spot size of the imaging beam. This small step size ensured there would be no significant change in the B-scan location within one heartbeat cycle, enabling us to apply a moving window in the post-acquisition alignment, which maximized the efficiency of frame usage and the reconstruction quality (described in Section 3.5). As an example, we acquired 30,000 B-scans at 100 Hz over 1025 µm, leading to a linear step size of ∼0.034 µm. With a ∼2 Hz heart rate, the total imaging beam movement over one heartbeat cycle was ∼1.7 µm, smaller than the spot size of the imaging beam. Acquiring this dataset took 5 minutes, during which, embryonic development occurred at a small spatial scale which was not expected to affect the 4D reconstruction. In the event of a whole-sample movement caused by the environment during data acquisition, the dataset would be discarded. The structure images were obtained from the intensity of the A-scans, with the DC component removed, and were reported in decibels. Within each B-scan, Doppler images were computed using the relative phase of A-scans through the 10 × 10 windowed Kasai autocorrelation [49], and the axial blood flow velocity was color-coded with cyan/magenta corresponding to the axial direction, Y, with a linear brightness representing the flow speed along that axis.

 

Fig. 2. An illustration of the entire 4D post-acquisition synchronization process. (A-L) represent the critical steps cited in the texts when describing the synchronization method.

Download Full Size | PDF

3.3 Time-difference cardiogram (TDCG) for heartbeat profile

The 4D alignment algorithm centers on the concept of the time-difference cardiogram (TDCG), named to draw an analogy to the established process of gating with an electrocardiogram [50]. The TDCG is calculated as the 2D average of the pixel-based squared difference between adjacent B-scans (black line, Fig. 2(B)), computed as

Here, B is the pixel intensity at (x, y) of the B-scan image, while N and M are the B-scan width (along the X axis) and height (along the Y axis), respectively. Because the step size between adjacent B-scan locations is smaller than the spot size of the OCT imaging beam by about two orders of magnitude, the TDCG represents temporal changes of the beating heart rather than the changes induced by scan motion. Thus, the TDCG is a 1D profile of the change in B-scan intensity primarily over time, t, but secondarily over space, Z, and labeled as the B-scan index i. Obtaining this parameter is a simple process of averaging multi-dimensional data down into a 1D signal for subsequent synchronization, which dramatically lowers the requirements for the computer processor and memory.

Since the TDCG metric is fundamentally an average rate of change, many features in an OCT image of the mouse embryo can potentially contribute to it. It can be thought of as a superposition of the following phenomena weighted by the fraction of the imaging plane which they affect:

Many imaging planes tend to include multiple of these phenomena. Since the heart is the major source of motion at this stage of development, most of these phenomena are correlated, such as blood flow corresponding to the contraction of the heart. Additionally, since the heart tube of the mouse is curved, multiple cardiac structures usually intersect the same imaging plane, causing their movements to superimpose to produce one TDCG signal. This superposition produces a consistent relative spacing between landmarks, each of which is identifiable by its unique peak amplitude and profile, which depends on what expansion, contraction, or flow phenomenon it represents. We assume that the phase of the TDCG exactly corresponds to the cardiac phase, i.e., that there is negligible change in the apparent cardiac phase along the Z axis.

3.4 H-score for identifying B-scans containing the heart

The H-score is the product of three metrics (Fig. 2(C)), developed for automatically identifying the B-scans that contain the beating heart. All three metrics are functions of the B-scan index, i, which ranges from 1 to n, the total number of B-scans. The first metric, C, is defined as

This metric (yellow line, Fig. 2(C)) represents a simple polynomial curve designed to generate a relatively high weight for the center portion of the B-scans, based on the assumption that the heart is more likely positioned around the center of the 3D data acquisition rather than the edges. The second metric (blue line, Figs. 2(B) and 2(C)) represents the average amplitude of heart wall movements over the whole heartbeat cycle. It is obtained by low-pass filtering the TDCG and then rescaling to the range of 0-1. The third metric (purple line, Fig. 2(C)) characterizes the heart wall movement within the heartbeat cycle, shown as the variance of high-frequency components. A high-pass filter is first applied to the TDCG to remove the DC bias, followed by squaring to calculate the local variance and then low-passing the variance to yield the average trend of that local variance. This metric is similarly rescaled to the range of 0-1. The concept behind the second and third metrics is that the B-scans covering the heart contain relatively more cardiac movements over time, reflected from both the low-frequency and high-frequency components of the TDCG, and thus, higher values of these metrics indicate the cardiac region.

With these three metrics, the H-score (black line, Fig. 2(C)) is calculated as their product (representing a logical AND), which takes all three factors described above into account when estimating the range of B-scan indices of the TDCG corresponding to the beating heart. The span of B-scans containing the heart starts at the first H-score value exceeding an arbitrary threshold and ends at the last H-score value exceeding that same threshold (red dashed line, Fig. 2(C)). We arbitrarily set the threshold as the midpoint of the minimal and maximal H-score for each given heart.

3.5 Synchronization process for aligning the heartbeat phase

The synchronization process converts the input dataset that has a combined Z and t axis into an output dataset in 3D with heartbeat phase φ. This conversion is made possible by assuming the semi-periodicity of the heartbeat, meaning that for any phase of the cardiac cycle, the heart must return to that phase after a variable period of time. Because the embryonic heart beats under its natural condition and there is finite sampling in time, the period of the heartbeat, shown as the number of B-scans, varies slightly between beats. However, the 4D output needs to have a fixed size for each axis. Thus, the average period of heartbeat is first obtained in terms of the number of B-scans and is set as the output period (OP) of the 4D reconstruction. Autocorrelation of the portion of TDCG covering the heart, as identified with the H-score, is performed to determine the OP (Fig. 2(D)), which represents the length of one cycle of the TDCG for alignment.

The starting point of alignment is then selected, which is at the maximal TDCG value within the cycle of the TDCG that contains the maximal H-score (Fig. 2(E)). This assigns the fastest moving phase of the heart as the starting phase of the final 4D reconstruction and ensures that the initial cycle of TDCG is taken from B-scans that contain clear cardiac activity. This TDCG cycle is utilized for an initial alignment of its adjacent four cycles of the TDCG, two each from immediately lesser and greater Z coordinates (Fig. 2(F)), to obtain five aligned TDCG sequences (Fig. 2(G)), whose average sequence serves as the reference sequence for subsequent forward and backward local alignments to synchronize the entire dataset. Such initial and subsequent local alignments share the same principle of a 1D cross-correlation, from which the Z offset producing the maximal value of correlation is recorded. The moving range for the cross-correlation is set as the size of one OP plus/minus 10% of it, rounded up. Prior to locating the maximum, a low-pass filter is applied to the resulted correlation profile, and this is for disambiguating the selection of the maximal correlation value when a flat region over the maximum is present.

With the initial alignment producing the reference sequence, the subsequent local alignment moves in both Z directions to avoid directional biases (Fig. 2(H)). A running average of five aligned and adjacent TDCG sequences is performed to update the reference sequence, where a newly aligned sequence of TDCG is added and the oldest or furthest sequence is removed (Fig. 2(I)). This improves the continuity between adjacent cycles by aligning them to an average of their neighbors; however, when the neighboring cycles correlate poorly with the reference sequence, adding them distorts the reference sequence. Thus, the adding of a new cycle is only conducted when the maximum of its cross-correlation with the reference sequence exceeds a threshold that is defined arbitrarily as 0.6 in our study. In the case where the maximum correlation is below the threshold, the Z offset is simply recorded as one OP, essentially retaining the previous alignment phase. This is necessary for TDCG cycles containing minimal or no movements induced by the heartbeat, usually the region relatively far from the heart. The use of a running average for the reference sequence with a control based on the correlation prevents the accumulation of phase drift in the Z direction. This process is summarized in Fig. 2(J).

With the Z offset determined for all TDCG cycles over the entire dataset, the φ coordinate is assigned to each sequence, with a range from 0 to OP-1 that corresponds to 0-100% of cardiac cycle. The B-scans are then written back to the computer drive and are labeled with the aligned Z and φ coordinates (Fig. 2(K)). Additionally, the same coordinates can be applied to any other image channels, such as Doppler B-scan images in our study (Fig. 2(L)). The whole algorithm is implemented in MATLAB (MathWorks), and the code is publicly available on our GitHub [51].

3.6 Experimental demonstration of the method

The scans used for developing and optimizing the 4D synchronization method were different from the ones used for testing and demonstration. To show the effectiveness and robustness of this method, we compared the result with an existing approach for the embryonic mouse heart [21,26,44] and tested the method with four mouse embryos featuring heart rates ranging from 1-3 Hz. Enabled by the high efficiency of this method, we reconstructed 4D cardiodynamics and hemodynamics over a 2-hour period of embryonic development, showing the growth and functional change of the mouse heart at a high spatiotemporal resolution.

3. Results

We first show typical intermediate data of the synchronization process (Fig. 3). This specific scan contains 30,000 B-scans that cover the entire beating embryonic heart. Figure 3(A) shows the TDCG across all 30,000 B-scans, and the blue line is the low-pass filtered signal representing the average movement amplitude over the whole heartbeat cycle. Zooming into the red boxed region where the maximum of the average amplitude is located, Fig. 3(B) shows a detailed view of the TDCG over four heartbeat cycles. The four windows are aligned to start with the maximal TDCG value, representing the starting phase for the alignment. Figure 3(C) shows the three metrics whose product provides the H-score used to estimate the ranges of B-scans that are likely to include the heart (Fig. 3(D)). Figure 3(E) shows the autocorrelation of the portion of TDCG above the threshold in Fig. 3(D) used to determine the OP (measured in the number of B-scans), which is 48 for this example scan. Figure 3(F) shows the maximum cross-correlation from each local alignment as a function of the Z coordinate in the final 4D output. The spike around the center with a correlation of 1 belongs to the initial TDCG cycle for starting the alignment, and the red dashed threshold line indicates the 0.6 cutoff for a correlation to be considered valid.

 

Fig. 3. Typical intermediate data of the synchronization method showing the core of the method and its self-monitoring parameters. (A) TDCG and its low-pass filtered signal. (B) Zoomed-in view of the TDCG over four cycles from the red boxed region in (A). (C) Three metrics produce the H-score. (D) H-score with the threshold determining the location of the heart. (E) Autocorrelation to estimate the OP. (F) The maximum of cross-correlation for alignment. (G) The difference in the original Z indices of B-scans that are equal in phase in the final output, and (H) its histogram. (I) The physical change in location of B-scans due to the alignment assigning a new Z coordinate, i.e., the cumulative distortion of the Z axis.

Download Full Size | PDF

The last three panels in Fig. 3 (G-I) constitute metrics of the synchronization output, used for insight into the alignment process. One critical parameter (Fig. 3(G)) is the difference in the original Z indices of B-scans that correspond to the same phase in the output image. This indicates the actual period of the heartbeat as perceived by the alignment process. As shown in Fig. 3(G), the values appear to be within 48 ± 1 B-scans, presenting a small deviation from the OP of 48 estimated based on the autocorrelation (Fig. 3(E)). The flat region to the left in Fig. 3(G) corresponds to the region in Fig. 3(F) where the correlation is less than the threshold, causing the Z offset to be directly assigned as one OP. Figure 3(H) shows the histogram of this parameter. Because the image is acquired as a linear 3D scan, each original B-scan has a distinct physical location where it was acquired. After the synchronization, B-scans within one aligned sequence are treated as if they are from the same physical location and are assigned a new Z coordinate, introducing distortion in the Z direction. Thus, we measure another parameter: the physical distance between the new Z coordinate of an aligned sequence and the original Z coordinate of the first B-scan in that sequence (Fig. 3(I)). This characterizes the cumulative physical distortion of the alignment, which, in this example, is smaller than one pixel and is well below the spot size of the imaging beam (∼4 µm), indicating that the distortion is negligible.

The comparison of the new synchronization method with the existing approach for the same dataset of mouse embryonic heart is shown in Fig. 4 and Visualization 1. The unaligned data serve as the reference to show the necessity of synchronization. The new TDCG-based method has a higher efficiency in utilizing the original B-scans, leading to a better spatial sampling in Z direction in the final 4D output. Thus, an overall improved imaging quality can be seen from the TDCG-based method (Figs. 4(A) and 4(B), 4(A’) and 4(B’)). In terms of the quality of the alignment, the TDCG-based method shows smoother 3D structures without distortions in comparison with the existing approach. From Fig. 4, these are specifically shown in the yolk sac (Box 1), the myocardium (Boxes 2 and 3), and the endocardium (Box 3). These indicate that the TDCG provides a superior sensitivity in identifying the heartbeat phase and that the new alignment method reliably synchronizes the phases from B-scans at different locations in the heart.

 

Fig. 4. Comparative assessment of the new TDCG-based method with the existing approach for mouse embryonic heart shows a superior reconstruction quality from the TDCG-based method (see Visualization 1). The 3D snapshots are shown at two consistent time points for each method. (A and A’) the TDCG-based method and (B and B’) the existing approach. (C and C’) The unaligned data are shown as a reference. Boxes 1-3 highlight the differences, and the arrow in C’ points at the significant distortion without alignment. The scale bars are 100 µm.

Download Full Size | PDF

The robustness of the method for synchronizing mouse embryonic hearts at different stages of development and, particularly, of different heart rates, is demonstrated in Fig. 5 and Visualization 2, Visualization 3, Visualization 4 and Visualization 5. Our data include embryos at four stages of cardiac looping, featuring heart rates of 2.13 Hz, 1.74 Hz, 1.52 Hz, and 1.24 Hz. En-face fly-through visualizations of the 4D cardiodynamics show the well-synchronized heartbeats, indicating the method can be used to image the early mouse heart at various developmental stages during dextral looping. Notably, at the stage of early E8.5, when the heart just starts to have coordinative beats with only few blood cells in circulation and relatively thin cardiac structures, the method can still successfully reconstruct the 4D heartbeat, clearly showing the time delay of contractions from the atrium to the ventricle (Fig. 5(D)). Because multiple cardiac cycles are used to reconstruct one heartbeat in the post-acquisition synchronization, the number of blood cells shown in the 4D visualizations (Fig. 5) appear to be greater than in reality.

 

Fig. 5. Imaging and reconstructions of 4D cardiodynamics from mouse embryos with different heart rates: (A) 2.13 Hz, (B) 1.74 Hz, (C) 1.52 Hz, and (D) 1.24 Hz. (See Visualization 2, Visualization 3, Visualization 4 and Visualization 5 for the hearts in A-D, respectively). The scale bars are 100 µm. A: atrium, V: ventricle, IF: inflow, LV: left ventricle, RV: right ventricle, DA: dorsal aorta.

Download Full Size | PDF

From the same hearts presented in Fig. 5, we show the integration of Doppler contrast with the 4D cardiodynamic images, where quantitative intracardial blood flow information is revealed in 4D together with the structural information of the heart (Fig. 6 and Visualization 6, Visualization 7, Visualization 8 and Visualization 9). Snapshots in Fig. 6 at specific time points show forward flows in the atrium-to-ventricle direction. Apart from the early E8.5 heart in Visualization 9, all hearts (Visualization 6, Visualization 7 and Visualization 8) experience retrograde flows located at the end of systole and/or at the beginning of diastole. As the retrograde flow is an important factor in the heart valve development, combining structural and functional information in 4D synchronization provides the opportunity to investigate how particular flow patterns are produced by heart wall dynamics and creates dynamic phenotyping to potentially study the role of flow-induced shear force in various aspects of early cardiogenesis.

 

Fig. 6. Combined 4D cardiodynamic and hemodynamic reconstructions of the mouse embryonic hearts beating at different rates: (A) 2.13 Hz, (B) 1.74 Hz, (C) 1.52 Hz, and (D) 1.24 Hz. (See Visualization 6, Visualization 7, Visualization 8 and Visualization 9 for the hearts in A-D, respectively). The scale bars are 100 µm. A: atrium, V: ventricle, IF: inflow, LV: left ventricle, RV: right ventricle, DA: dorsal aorta, OFT: outflow tract.

Download Full Size | PDF

In Table 1, we show the statistics of the synchronization processes performed on the four embryonic hearts presented in Fig. 6, highlighting the efficiency of this new method. The computer running the algorithm has a Core i9-10900X CPU and DDR4 2933 MHz RAM, while the reading and writing of files is based on a typical solid-state drive through the SATA III interface. All data is uncompressed, and the B-scan utilization shows what fraction of the original B-scans are used in the aligned 4D output. The ∼100% utilization of B-scans is achieved due to the design of the method that allows for grouping any continuous portion of the original B-scan series as a full cardiac cycle for alignment. Greater than 100% utilization (such as 101.1% from Embryo B) is possible since the aligned TDCG sequences can overlap, meaning the aligned time sequences from adjacent Z locations of the 4D output could share the same original B-scans at the beginning or the end of the sequences. The total data processing time is broken down into the file reading time, alignment time, and file writing time. Here, file reading includes fetching images out of the RAW file of structure OCT data and simultaneously converting them to the TDCG to reduce memory overhead. Alignment starts with estimating OP and ends when all the Z offsets are determined. File writing time involves reading structure and Doppler OCT images from the corresponding RAW files one by one and rewriting them as individual TIFF files with their corresponding final 4D Z and φ coordinates in name. Because the alignment process only relies on the 1D TDCG profile, sub-second alignment time is achieved with a mid-level CPU, while the file reading and writing time is primarily determined by the hard drive interface.

Table 1. Statistics of the synchronization process

View Table

The efficiency of this method allows us to reconstruct many datasets within a reasonable time, enabling repeated imaging and visualization of the mouse embryonic heart development with 4D cardiodynamics and hemodynamics (Fig. 7 and Visualization 10, Visualization 11, Visualization 12 and Visualization 13). Visualization 10 and Visualization 11 show the cardiovascular development over two hours starting around E9.0, with the cross-sectional views on the atrioventricular and interventricular regions, respectively. The two-hour developmental time was sampled at a 12-minute interval. The growth of the heart, the vein, and the aorta, as well as the change in the blood flow, are directly revealed from the visualizations. These can also be seen from the starting and ending time points shown in Fig. 7. The duration and the cycle fraction of the retrograde flows at the inflow and outflow regions over time are respectively shown in Figs. 7(C) and 7(D), where a decreasing trend at the outflow is clear, likely contributed by a gradual change of the ventricular wall dynamics as the heart remodels. From Visualization 12 and Visualization 13, the cardiac looping is already significant at the start of imaging, while the expansions of the atrium and ventricles over two hours produce further bending of the heart. Such repeated 4D imaging over development could bring new insights into the dynamics of cardiogenesis at extended time scales.

 

Fig. 7. Dual-contrast 4D imaging and reconstruction of the mouse embryonic heart development over two hours. Snapshots at (A and A’) the atrioventricular region and (B and B’) the interventricular region at the first and last developmental time points. Plots of (C) the duration of the retrograde flow and (D) its ratio to the heartbeat cycle at the inflow and outflow regions over development. The scale bars are 100 µm. A: atrium, RA: retrograde flow, IF: inflow, LV: left ventricle, RV: right ventricle, OFT: outflow tract.

Download Full Size | PDF

4. Discussion

Previous post-acquisition synchronization methods use B-M scans [35–38], which generally require more B-scans than needed to select data subsets for alignment of the heartbeat phase. In comparison, the new method utilizes a densely sampled, linear 3D scan; it not only simplifies the data acquisition process that is available from most commercial OCT systems but also leads to a ∼100% B-scan usage. Notably, this allows the aligned sequences to overlap when the heart rate in the corresponding cycle is slightly lower than the average (OP). Relying on this scanning scheme produces two adverse effects related to the spatial scale. One is a small amount of scale distortion in the Z direction of the final 4D output, which is characterized and shown in Fig. 3(I). This results from treating a nonlinear Z coordinate as a linear one. Although it is possible to correct this by resampling the aligned 4D onto a linear Z coordinate, in practice, the distortion is such a small portion of a pixel for each B-scan that even their accumulation is far less than one pixel (and the spot size of the imaging beam), indicating that a correction is unnecessary. The other adverse effect is a shift of the aligned 4D in the Z direction over time, cumulating in exactly one pixel of offset, which is not visible during the heartbeat but can be noticed as a flicker when the 4D visualizations wrap back around from the end to the beginning. This Z-direction shift comes from the scanning (movement of the B-scan location) over one heartbeat cycle during data acquisition. Since the total distance traveled during a heartbeat is smaller than the transverse resolution of imaging, the total distance of the one-pixel shift is also smaller. The shift could potentially be further reduced by having an even smaller sampling interval along the Z direction of 3D acquisition, or it could be removed by resampling the frames using interpolation to eliminate the fraction of a pixel shift in each aligned sequence.

The synchronization is based completely on the pattern of cardiac movement. Thus, when the TDCG cycle containing B-scans that do not cover any portion of the beating heart, the local alignment likely produces low correlation values, typically shown as the region below the threshold line in Fig. 3(F). In such regions, the OP is assigned as the Z offset, and therefore, the blood flow in vessels can appear poorly synchronized. Since blood cells largely do not return to the same locations at the same phase of the heartbeat, structure-image-based correlation might not have sufficient sensitivity to construct a heartbeat profile based on minor tissue movements induced remotely by the beating heart. The blood flow velocity information from the Doppler OCT data could potentially be used to assist the alignment in these regions.

Since this alignment method is based on changes in the image over time, it can yield incorrect phase alignments when movement patterns between B-scans are identical except for a phase offset in time, as is the case when the imaging plane is exactly perpendicular to the axis of a peristaltic tube. Fortunately, this issue can be avoided by orienting the B-scan plane to intersect the longitudinal cross-section of the tube instead, such that the phase of the contraction wave is almost constant as the plane moves along its normal. In the case of the embryonic mouse heart during dextral looping, the issue is absent, since the tube’s curvature forms a spiral, and a transverse section will always intersect the imaging plane a second time. When multiple parts of the cardiac tube or other structures intersect the B-scan plane, the opposing direction of the two peristaltic waves prevent the apparent cardiac phase from changing as the B-scan plane progresses along its normal. Additionally, the large fraction of each B-scan occupied by the mouse myocardium and other tissue peripheral to the heart generates a low-frequency baseline TDCG signal which is relatively constant in phase throughout the entire acquisition. Our correlation-based approach allows us to consider multiple simultaneous movement phenomena to enhance our phase reliability beyond that of simple maximum peak detection.

We usually calculate the TDCG over the whole B-scan image. Although this reduces the TDCG signal-to-noise ratio (SNR) by including regions of significant background, which contain high speckle variance while lacking visible cardiac structures, the SNR is sufficient for alignment of the dataset whenever the heart occupies a significant portion of the B-scans. Otherwise, we can improve the SNR by specifying a rectangular subset shared across B-scans (a crop of the image) for calculating the TDCG. We only use this crop to calculate the TDCG; we do not retain it in the aligned image.

As shown in Table 1, the time used for alignment is negligible compared with the time needed to read and write images. Upgrading the data storage interface from SATA to NVMe would accelerate the file reading and writing processes, but we could further reduce the IO time by adjusting the strategy of handling the data and visualization. File read speed could be improved by multithreading the TDCG computation step to compute the previous difference of images while the following image loads in parallel. On the output side, the third-party visualization software (Imaris) requires us to rewrite images as individual files with Z and φ coordinates in the names, which is the most time-consuming step. Building and applying a customized visualization software that directly reads images from the original RAW file in the order specified by the Z offsets would eliminate that step.

Our algorithm currently requires the user to enter an initial estimate of the heartbeat period. It uses this estimate to disambiguate between integer multiples of the OP obtained during autocorrelation. E.g., when the true OP is 12 B-scans, the periods 24 and 36 should correlate equally well, in theory. This estimated heart rate also guides the specification of filters that generate two of the three metrics for calculating the H-score. Potentially, this user-defined step can be omitted by performing autocorrelation over a wide range of possible cardiac periods and employing additional signal processing techniques to disambiguate integer multiples of the smallest period correlated, which can further improve the automation of the method.

In this work, our B-scan acquisition rate is either 100 Hz or ∼83.3 Hz, corresponding to the volume rate of the final 4D visualizations. For the mouse embryonic heart, this imaging speed is sufficient and favorable for studying the temporal characteristics of the heartbeat. We designed the low-pass filter to smooth the correlation during local alignment over this specific range of B-scan rates. Specifically, the filter taps are [0.1875, 0.6250, 0.1875], corresponding to a normalized cutoff frequency of approximately 0.43 B-scan^-1. When using substantially different B-scan acquisition rates, it may be necessary to adjust the taps of this filter.

An important feature of our method is that the starting heartbeat phase of the final aligned 4D is not random but is the phase with the largest (or fastest) heart wall movement. Phase stability is beneficial for reconstructing the repeated 4D imaging of the heart development (Visualization 10 and Visualization 11) since the 4Ds are automatically connected approximately in phase, making potential quantitative analyses convenient. However, the answer to which phase of the heart cycle corresponds to the fastest motion (maximal TDCG) is unstable over extended developmental periods. Therefore, a further, more precise alignment of heartbeat phases among the 4Ds should ideally be conducted based on another layer of TDCG assessment.

In summary, we report a new post-acquisition synchronization method for 4D OCT imaging of the beating mouse embryonic heart. This open-source method uses a single, dense, linear 3D scan of the sort readily available in most OCT systems. We demonstrated its 4D reconstruction quality and robustness for a wide range of heart rates by comparing an existing approach and demonstrating synchronization of scans taken at different developmental stages. The method features a ∼100% B-scan usage and sub-second alignment time. Enabled by this high efficiency, we show repeated 4D imaging and visualizations of the mouse embryonic heart development over two hours. This method could enable a broader application of OCT for dynamic, 4D analysis of the early cardiogenesis.