Semi-automated shear stress measurements in developing embryonic hearts

: Blood-induced shear stress inﬂuences gene expression. Abnormal shear stress patterns on the endocardium of the early-stage heart tube can lead to congenital heart defects. To have a better understanding of these mechanisms, it is essential to include shear stress measurements in longitudinal cohort studies of cardiac development. Previously reported approaches are computationally expensive and nonpractical when assessing many animals. Here, we introduce a new approach to estimate shear stress that does not rely on recording 4D image sets and extensive post processing. Our method uses two adjacent optical coherence tomography frames (B-scans) where lumen geometry and ﬂow direction are determined from the structural data and the velocity is measured from the Doppler OCT signal. We validated our shear stress estimate by ﬂow phantom experiments and applied it to live quail embryo hearts where observed shear stress patterns were similar to previous studies.


Introduction
Blood flow is a critical epigenetic factor during the early stages of cardiac development. Altered hemodynamics are thought to contribute to the phenotypes associated with congenital heart defects (CHDs) which affect ∼1% of births worldwide [1]. Several studies have shown that disturbing blood flow patterns via surgical or stimulation techniques results in malformations to key cardiac structures which leads to CHDs [2][3][4][5][6][7]. Altered blood flow patterns create abnormal shear stresses on the endocardium of the heart tube. These stresses are transduced which activates intracellular signaling pathways that lead to alterations in gene expression [8,9]. Over the last several years, researchers have demonstrated that altered blood flow leads to changes in gene expression in looping embryonic hearts [10][11][12]. We are beginning to understand how blood flow impacts early heart development, but better tools for measuring and analyzing shear stress are needed.
Three optical coherence tomography (OCT) approaches have been developed for evaluating shear stress in the early-stage embryonic heart. The first approach uses a combination of modeling and imaging. Liu et al., [13,14] developed image-based finite element (FE) models of the outflow tract (OFT) in a chick embryo where 4D OCT images were used to extract the geometry of the wall over time and blood pressure was measured using a transducer at the inlet and outlet of the OFT to provide the boundary conditions. In another method by Midgett et al., [15], velocity was computed from longitudinal B-mode images assuming a parabolic velocity profile, which only enabled calculation of the average wall shear stress. The other approach measures the shear stress directly from Doppler OCT blood velocity measurements. Peterson et al., [16] mapped the shear stress in hearts of quail embryos directly using 4D OCT images where the wall geometry was obtained from structural images and the blood velocity was measured from the Doppler signal.
These approaches are exciting and have great potential, but they involve complex instrumentation, intensive processing or are not able to provide the actual variations of wall shear stress on the endocardium. Simpler methods to assess shear stress would make this technology useful to more researchers and increase its utility for biological studies.
Longitudinal cohort studies of cardiac development involve monitoring a large number of embryos over the span of several weeks, requiring fast data throughput as well as user-friendly processing methods. As a result, a more practical approach is needed to incorporate shear stress measurements in longitudinal studies. In this manuscript, we introduce a new method to quickly estimate the shear stress of a selected OCT frame without the need to capture volumetric data. Coordinates of an adjacent frame help determine the direction of flow measured by Doppler OCT and the angles needed to project the OCT image to an orthogonal plane of the heart tube to more accurately assess shear stress. We validated this shear stress measurement method by flow phantom experiments, and applied it to live quail embryo hearts.

Method of estimating shear stress
Shear stress τ is calculated using the equation τ = ηdv/dr where η is the viscosity, v is the velocity measured from the Doppler OCT signal and r is the distance from the wall along the normal direction. A simplified sketch of the tubular heart imaged by an OCT system is shown in Fig. 1(a). The OCT image is not perpendicular to the direction of blood flow in the heart. Shear stress requires measuring the velocity gradient normal to the wall (dv/dr), but the image does not contain this information. However, because flow is laminar in the embryonic heart causing shear stresses to vary slowly on small spatial scales, one can assume that the velocity gradient can be calculated by finding a new distance dr that corresponds to dr of the closest orthogonal plane and allows one to find v in the OCT image and approximate the velocity gradient using v and dr ( Fig. 1(b), zoomed in box). The location of the points on the orthogonal cross section can be found by projecting the luminal wall in the OCT plane onto the orthogonal plane. The detailed description of our shear stress measurement is given in section 3.3. Fig. 1. Simplified sketch of a tubular heart imaged by an OCT system: a) the side view of the heart where blood, inner and outter curvatures are indicated, b) the image plane is parallel to the OCT beam and shear stress is measured on the orthogonal cross section. Zoomed box: the point at a distance dr and velocity of v on the orthogonal cross section corresponds to a point located at a distance dr on the imaging plane. A distance dr on the imaging plane would give an incorrect velocityv oct .

Embryo preparation
IACUC approval was not required for this study which involves the use of quail embryos that were collected at embryonic day 2 at the latest. Fertilized quail eggs were incubated in a humidified incubator (Eppendorf New Brunswick, Germany) at 38°C. At 48 hours of development, the eggshells were removed and the embryos were cultured in Petri dishes. Tubular hearts of the embryos were imaged at 48 hours of development.

OCT imaging
The high-speed OCT system used for data acquisition includes a Fourier Domain Mode Locked swept laser source (Optores GmbH, Germany) operating at a rate of 1.68 MHz as previously described [17]. The OCT interferometer was built in a Mach-Zehnder configuration, where the sample arm was placed inside an incubator to control the temperature and humidity during imaging of live quail embryos. The axial resolution is ∼12 µm in air and the beam spot size at the imaging focal plane is approximately 15 µm over a 4×4 mm 2 field of view. Adjacent B-scans were acquired over a lateral length of 1 mm at a step of 10 µm. Each B-scan consisted of 16000 A-lines with 594 pixels in depth (acquisition time 9.52 ms). OCT images were collected over multiple heartbeats at sequential slice locations and reordered using image-based retrospective gating [18].

Shear stress measurement
Shear stress τ is calculated using the equation τ = ηdv/dr where η is the viscosity, v is the velocity and r is the distance from the wall along the normal direction. To measure shear stress in the tubular heart, we made similar assumptions as previous papers [16]. Assumptions include blood being a Newtonian fluid with a dynamic viscosity of 5 mPa s, blood flow in the tubular heart having low Reynolds and Womersley numbers indicating laminar flow dominated by viscous forces, and blood flowing in the direction of the center line of the heart tube. Because the blood flow is laminar, we also assumed that shear stresses vary slowly across small spatial scales. Our method of estimating the shear stress includes the following steps which were implemented in Matlab R2016a (MathWorks Inc., USA): 1. Segmenting the lumen and determining normal directions to the lumen: To determine the location of the wall on which the shear stress is calculated, the lumen is segmented manually from the OCT structural image. We assume the OCT image is in the x-z plane, where the z axis corresponds to the axial depth (along an A-scan) and x corresponds to lateral dimension cutting across the A-scans. Using the roipoly function which creates an interactive polygon selection tool, a few points are selected on the luminal wall and a binary image is returned with the pixels inside the selection region set to 1. The points on the boundary of the selected region are found by the boundary function. To define the luminal wall (l), the cscvn function fits a spline curve through the boundary points and a smoothing filter is applied by the sgolayfilt function. Having defined the luminal wall, the normal directions are calculated. The slope of the normal line is the negative reciprocal of the slope of its tangent at any point around the curve. The derivatives of l, the curve describing the luminal wall with respect to x (dl x ) and z (dl z ) are found by the diff function to define the slope m of the normal line: 2. Finding flow direction and orientation of the orthogonal cross section: In order to find the direction of flow, an adjacent frame to the selected frame is needed. The lumen centroid of each frame is found by using the mean function on the luminal wall points. By finding the line that connects these two centroids, one can identify the flow direction. The cross section on which the shear stress is measured is defined by the plane orthogonal to the flow vector. Considering the OCT image in the x-z plane, the orientation of the plane orthogonal to the flow direction is given by a polar angle α p and an azimuthal angle α a (Fig. 2(a)).
3. Projecting the OCT plane onto the plane orthogonal to the heart tube: Each point l i on the luminal wall in the OCT plane is projected onto the orthogonal plane to find the corresponding point w i . The projection of the points in three dimensions is achieved by two sequential two-dimensional projections using angles α p and α a (Fig. 2(b) and (c)). Figure 2(b) shows the projection in the y-z plane using α p only. The projected z components are defined by: This projection results in a narrowing of the OCT image along the z direction while the x components remain unchanged. Figure 2(c) shows the projection in the x-y plane which uses α a only. The projected x components are defined by: This projection results in a narrowing of the OCT image along the x direction while the z components remain unchanged.
4. Determining the position of l i : The point w i located at a selected distance dr from w i corresponds to l i which is located at distance dr i along the normal direction to the luminal wall. In the OCT plane ( Fig. 2(d)), each pair of w i and l i is defined by a line where the slope is the projected m i (slope of the calculated normal at w i ) : The distance between each pair of w i and l i is the constant dr: Using Eq. (3), we can solve Eq. (4) for x l i and z l i : Although the selected dr is constant for all the points on the orthogonal cross section, the distance dr i on the OCT plane varies from point to point around the lumen. The example points l 1 -l 3 in Fig. 2(d) illustrate the change in dr i from the maximum value dr 1 to the minimum value dr 3 .

5.
Measuring the flow velocity using complex regression Doppler: A thorough shear stress measurement relies on accurate quantification of flow velocity in the vicinity of the wall. Since blood flow and wall motion undergo significant increases in velocity during heart development, imaging parameters need to be adjusted continuously to accommodate the large velocity range required in longitudinal cohort studies of embryos. Doppler phase shift methods that enable detection of slow velocities, result in phase wrapping of faster flows at the center of the tube. We used complex regression Doppler OCT. The recently published method for measuring Doppler phase shift extends the dynamic range by lowering the minimum detectable velocity while providing higher accuracy and precision [17]. 6. Calculating shear stress: The velocity v i at point l i corresponds to the velocity at a distance dr on the orthogonal cross section or dr on the OCT plane. The v i is measured from the Doppler image and is corrected by the Doppler angle (α p ) to find the absolute velocity v abs . Shear stress is estimated by dividing the change in absolute velocity dv abs over the distance dr and multiplying that by viscosity η.

Phantom validation experiment
To evaluate our method of estimating shear stress, 2% lipid solution (Intralipid, Clayton, NC) was pumped through a capillary tube with an inner diameter of 300 µm using an NE-300 Just Infusion syringe pump (New Era Pump Systems Inc., USA). The flow rate was set to 25 ml/h and the tube was immersed in water to avoid aberrations. The images were acquired at varying Doppler angles (55°to 80°). The angles were confirmed from structural OCT volumes at each position.

Results
Flow phantom experiments were performed to evaluate the performance of our shear stress estimate. The structural image of a cross section is shown in Fig. 3(a) where the selected points on the boundary (blue dots) were used to define the luminal wall (white curve). The centroid of this luminal wall and that of the adjacent frame were used to determine the flow direction which was very similar to the direction found using several frames. A Doppler OCT image can be seen Fig. 3(b) where the luminal wall is indicated by a white curve. The black and dotted black curves correspond to the points located at dr = 45 µm and corrected dr from the wall, respectively. The numerical variation of dr (red points) around the lumen is plotted in Fig. 3(c) where dr (blue points) remains constant. The comparison of measured shear stress using the velocities at dr and dr can be seen in Fig. 3(d) where the black dotted line indicates the expected value of 2.2 Pa.  Figure 4 shows the shear stress along a capillary tube positioned at an angle of 63°. When using measured velocity at dr (uncorrected) for dr = 45 µm and η = 1 mPa s [16], the shear stress at the top and bottom of the tube appears lower than the side surfaces where the expected value of 2.2 Pa is observed (Fig. 4(a)). Measured velocity at dr (corrected) results in a uniform shear stress around the cross sections ( Fig. 4(b)). To further validate our method, the Doppler angle was varied from 55°to 80°and the shear stress was measured at dr = 45 µm and the corresponding dr . At each angle, the shear stress was measured along the periphery of the capillary tube on five frames and all the measured values were averaged. Figure 5 illustrates the mean of the measured uncorrected shear stress (blue bars) and the mean of corrected shear stress (red bars) along with the standard error of mean for the sample size of 75. While uncorrected shear stress is variable with larger standard error, corrected shear stress is consistently close to the expected shears stress of 1.05 Pa and also has a smaller standard error. As the Doppler angle increases, the OCT plane becomes closer to the orthogonal cross section, however, the SNR of the Doppler signal becomes weaker, which may explain the lack of an apparent trend between the two measurements. An example of estimating the shear stress in a day 2 quail embryo is shown in Fig. 6. A high quality, averaged image along the tubular heart is illustrated in Fig. 6 (a) to provide context for the location of the following images where the inner and outter curvatures are indicated and represent the approximate location of the sebsequent cross section. The structural image of a cross section overlaid by the Doppler signal is shown in Fig. 6(b). The luminal wall was segmented by selecting a few points on the boundary (Fig. 6(c), blue dots) and velocities were measured at dr from the lumen (Fig. 6(d), black dotted line) on the Doppler OCT image. The estimated shear stress using dr = 45 µm and η = 5 mPa s [19] is shown in the inset of Fig. 6(b). Higher shear stress is observed on the inner curvature of the tubular heart which is consistent with the results reported previously from mapping the shear stress of the entire heart [16]. Here, we presented the wall shear stress at one time point, while it can be calculated over the entire cardiac cycle (shown in Visualization 1).

Conclusions
We introduced a new approach to estimate shear stress in an OCT frame by using only two adjacent frames, eliminating the need to process the entire volumetric data. This approach is quick, user friendly, and allows high data throughput. Flow phantom experiments confirmed the accuracy of the shear stress estimate. The previously published method of measuring shear stress in a developing heart [16] took several days whereas our method cuts down the processing time to a few minutes and results in similar shear stress patterns. This semi-automated approach can become fully automated by combining speckle variance OCT with gated imaging for segmentation of the luminal wall [20]. Although this method does not provide shear stress throughout the entire heart, well positioned B-scans can impart excellent insight into whether shear stresses are abnormal and affecting heart development. This practical method can be used to incorporate shear stress measurements in longitudinal cohort studies where a large number of embryos are monitored over several developmental stages.

Disclosures
The authors declare that there are no conflicts of interest related to this article.