Defining Coronary Flow Patterns: Comprehensive Automation of Transthoracic Doppler Coronary Blood Flow

The coronary microcirculation (CM) plays a critical role in the regulation of blood flow and nutrient exchange to support the viability of the heart. In many disease states, the CM becomes structurally and functionally impaired, and transthoracic Doppler echocardiography can be used as a non-invasive surrogate to assess CM disease. Analysis of Doppler echocardiography is prone to user bias and can be laborious, especially if additional parameters are collected. We hypothesized that we could develop a MATLAB algorithm to automatically analyze clinically-relevant and non-traditional parameters from murine PW Doppler coronary flow patterns that would reduce intra- and inter-operator bias, and analysis time. Our results show a significant reduction in intra- and inter-observer variability as well as a 30 fold decrease in analysis time with the automated program vs. manual analysis. Finally, we demonstrated good agreement between automated and manual analysis for clinically-relevant parameters under baseline and hyperemic conditions. Resulting coronary flow velocity reserve calculations were also found to be in good agreement. We present a MATLAB algorithm that is user friendly and robust in defining and measuring Doppler coronary flow pattern parameters for more efficient and potentially more insightful analysis assessed via Doppler echocardiography.

The coronary microcirculation (CM) is unique in physiologic structure and function compared to other microvascular beds. It plays a critical role in the regulation of blood flow and nutrient exchange to support the viability of the heart. Because this vasculature lies within the surrounding myocardium, it constantly experiences transmural forces in the form of systolic contraction and diastolic relaxation with every cardiac cycle. This dictates the amount of blood flow regulation to a large extent, in conjunction with the CM's response to neural and hormonal factors. Compared to other blood flow in the body, coronary flow is unique in that it primarily occurs during myocardial diastole, whereas the CM is mostly occluded during systole 1 .
Transthoracic Doppler Echocardiography (TTDE) has proven to be a useful and relatively inexpensive tool to non-invasively assess cardiac perfusion by measuring parameters such as coronary flow velocity reserve (CFVR) and velocity-time integral (VTI) [2][3][4][5][6] . For both human and animal subjects, blood flow is measured via one of the major coronary arteries (left main, right main, left anterior descending) at various windows under both baseline and stress conditions, which yields the characteristic biphasic coronary flow pattern (CFP). In many disease states however, the CM becomes structurally and functionally impaired leading to changes in TTDE measurements and can thus be indicative of progressing CM pathology.
For example, in adult humans with type II diabetes mellitus (T2DM), multiple studies have demonstrated a reduction in coronary flow reserve (CFR) or CFVR compared to matched controls despite being asymptomatic [7][8][9] . Previous studies by our laboratory have shown that the CM in young T2DM murine and porcine models undergoes early inward hypertrophic remodeling 10,11 , and others have demonstrated functional deficits 12,13 . Structural remodeling was associated in vivo with a reduction in coronary blood flow (CBF) and CFR. Moreover, this remodeling occurred prior to occlusive macrovascular atherosclerosis which further suggests the importance of examining the CM especially early in disease progression 11 . Of additional interest, recent studies have shown Transthoracic Doppler Echocardiography. Coronary blood flow velocity was measured noninvasively with a high-frequency, high-resolution ultrasound unit (Vevo2100, Visual Sonics, Toronto, Canada) equipped with a 30 MHz probe, at baseline (1% isoflurane), and under conditions of maximum flow (hyperemia, 3% isoflurane) as previously described 10 . Doppler measurements of the left main coronary artery diameter and flow were performed under a modified four chamber view. Mice were anesthetized with 2% isoflurane vaporized with 100% oxygen. Following induction, isoflurane was reduced to 1% to determine baseline coronary flow, and then increased to 3% to measure maximal coronary flow. Baseline and hyperemic PW Doppler files were either manually analyzed on the Vevo 2100 software or exported for automated analysis using the MATLAB program.
Coronary Flow Pattern Program. The CFP program was developed using MATLAB software (The MathWorks Inc., Natick, MA) and was designed to analyze and extract time intervals, velocity points, and slopes from PW Doppler CFP AVI video files exported from Vevo 2100 software. A detailed description of our algorithm for processing the video files and extracting the desired parameters can be found in the online supplemental section.
Data Exportation. Raw AVI video files were directly exported from the Vevo2100 software for analysis. Specifically, uncompressed PW Doppler AVI video files were automatically exported at a size of 880 × 666 pixels, a frame rate of 30 frames per second, and sweep speed parameter of 0.85 seconds.
Region of Interest Extraction. PW Doppler Files: The raw uncompressed PW Doppler video files contain the desired coronary Doppler window region as well as a time axis, a velocity axis, an ECG recording, a B-Mode window, and various study labels (Fig. 1A). In order to extract the complete coronary Doppler region and ECG recording, a number of cropping and parsing steps were initiated and were built upon similar techniques by Magagnin et al. 14 and Zholgharni et al. 20 . The final resulting images were the full Doppler CFP sequence from the zero-velocity baseline to the maximum velocity and the full ECG recording (Fig. 1B).
Feature Extraction and Parameter Measurements. PW Doppler Files -Envelope Overlay: After cropping the PW Doppler files to obtain both the full CFP sequence and the ECG recording, several image processing techniques were used to create an envelope overlay of the CFP. A flow chart of the image processing techniques is presented in Fig. 2.
PW Doppler Files -Parameter Extraction: The ECG recording was used to separate individual CFP cycles per heartbeat. Specifically, the start and end of each CFP cycle was identified as the R-peak to R-peak interval from the ECG. Following separation of each individual CFP cycle, the derivatives of their envelopes were calculated. Both the flow pattern envelope and its derivative were utilized in measuring our newly defined time intervals, velocity points, slopes, VTI and HR specified below (Fig. 3). The parameters were measured and stored for each cycle in the CFP sequence. In addition to extracting the specified parameters from PW Doppler files, the program averages each parameter per animal and calculates CFVR. Often in TTDE, certain cycles should be excluded from analysis, as they are not representative of cycles in the entire sequence ( Supplementary Fig. 1). Our MATLAB program also contains a function that eliminates non-representative cycles and outliers from analysis. Figure 4 displays a portion of a CFP analyzed with the program showing the identified velocity points and envelope overlay. Finally, CFVR was calculated by dividing average peak hyperemic velocity into peak baseline velocity: CFVR = PVhyperemia/PVbaseline. Validation Protocol. Two trained lab personnel with adequate experience in Doppler CFP analysis performed both manual and program analysis in a blinded manner. One PW Doppler baseline CFP file, and one PW Doppler hyperemic CFP file were analyzed per animal. 6 animals in total were analyzed. The PW Doppler files recordings were 1.20 seconds -4.90 seconds in duration with a minimum of 6 complete cardiac cycles.
Manual Analysis Protocol. Each rater analyzed the animal files a total of two times separated by two distinct viewings in a blinded manner. The individual PW Doppler files were blinded and randomized so the rater had no knowledge on the specific animal type or file order. The files were further randomized between viewings. Measurements were obtained in the Vevo 2100 software. Each rater analyzed the same number of CFPs as the start and end cycles were pre labeled prior to analysis. If the rater deemed any of the cycles to be unclear or non-representative, he/she noted that cycle as non-analyzable and continued analysis until the entire sequence was completed. The measurements obtained for each cycle included the 4 time intervals, 5 velocity points, 4 slopes, HR, and VTI. Once flow pattern cycle measurements were made in the Vevo 2100 software, the rater populated a pre-labeled Microsoft Excel datasheet. The total time taken to complete both viewings was recorded.
Program Analysis Protocol. The same raters performed the program analysis similarly in two separate viewings. The user ran the MATLAB program for each animal and the measured parameters and calculations were directly exported into a Microsoft Excel datasheet. The total time taken to complete both viewings was recorded as well as time taken to export the raw files into a folder for analysis.
Parameter Selection and Disease. After performing the validation study to determine if our program selected parameters similar to a trained rater, we utilized a larger sample size of mice to determine if any of our newly defined parameters differed between normal and disease states as a proof of concept. For this test, we analyzed CFPs with our program from 16-wk Db/db (n = 18) and db/db (n = 20) mice under baseline and stress conditions. The average and standard error of the mean (SEM) was found for each parameter and an unpaired t-test was used to determine if any significant differences were present between normal and diabetic CFP parameters.
Statistical/Data Analyses. Intra-observer variability, inter-observer variability, and variability between manual analysis and the MATLAB program were examined. For each of these tests, data was expressed as mean +/− SD. Bland-Altman analysis was executed to calculate bias (mean difference) with limits of agreement set as +/− 2 SD 21 . Linear regression analysis was performed to determine the coefficient of determination (R 2 ) and the regression equation. Unpaired t-tests with a significance level of p < 0.05 were also executed. Statistical tests were done in GraphPad Prism 7 software and Microsoft Excel. The CFP program was created in MATLAB and the program analysis protocol was run on a computer utilizing an Intel ® Core ™ i3-2100 CPU @ 3.10 GHz.

Results
For the PW Doppler files between all 6 mice, there were a total of 98 complete CFP cycles measured at 1% isoflurane (baseline) that could have undergone analysis. There were a total of 117 complete CFPs measured at 3% isoflurane (hyperemia) that could have undergone analysis. Because each rater analyzed the same files in two separate randomized viewings, each frame and CFP had the same chance of being analyzed twice. Supplementary Table 1 specifies the total number of complete CFP cycles that could have undergone analysis for the PW Doppler files separated by animal. PV, HR, and VTI results under baseline and hyperemic conditions are presented below as they are the most frequently measured in clinical practice.

Intra-Rater Variability.
To determine the intra-rater variability for manual analysis and for the program analysis, only CFPs that were analyzed in both viewings per rater were included. If a CFP was analyzed in only one viewing, or not analyzed in either viewing, the cycle was excluded. This was deemed appropriate as the twice analyzed cycles were likely the cleanest and most representative cycles. Supplementary Fig. 2 displays representative images of both the manual analysis and the program analysis over the same set of cycles to highlight the variability between viewings. Bland-Altman and linear regression analysis were performed on each rater's consistent cycles for PV, HR, and VTI. Table 1 shows the intra-rater variability of these parameters for both manual analysis and the program analysis. PV, VTI, and HR intra-rater variability was reduced in the program analysis for both raters as evidenced by a bias nearing closer to 0 and smaller limits of agreement in the Bland-Altman analysis. Additionally, the coefficient of determination neared closer to 1 in the linear regression analysis. HR variability in the program analysis was nullified completely as this section of the program is fully automated. The decreasing trend in variability with the program was found in both baseline and hyperemic CFPs for the listed parameters. Inter-Rater Variability. Supplementary Tablesnual and program inter-rater variability, only cycles that were analyzed in both viewings and that were consistent among each rater were included. As these CFPs were analyzed twice by both raters, they were deemed appropriate to incorporate into this analysis. Average values between the two viewings were calculated for each rater and the total number of cycles to determine the inter-rater variability is shown in Supplementary Table 2. Bland-Altman and linear regression analysis were performed on the rater's averages for their consistent cycles for PV, HR, and VTI. Table 2 displays the statistical figures for the inter-rater variability. Variability between raters was reduced when using the program to analyze CFPs for PV, HR, and VTI at 1% and 3% isoflurane. This is supported by the bias's tending closer to 0 and smaller limits of agreement in the Bland-Altman analysis as well as larger coefficients of determination in the linear regression analysis.

Manual Analysis vs. MATLAB Program Analysis. While the MATLAB program reduced variability
individually and between different raters, it was equally important that the program selected parameters that were in good agreement with manual analysis. The program was tested to assess its ability in finding total parameter averages per animal.
Animal by Animal Validation. Often, parameters are averaged across a number of flow patterns for each subject/ animal in order to obtain representative measurements. The purpose of this validation was to compare the average parameter output of the program to the manual average for each mouse. For each CFP per animal at 1% and 3% isoflurane, the mean PV, VTI, and HR were calculated by averaging both raters' measurements. If a cycle was not analyzed by either rater, that cycle was excluded (Supplementary Table 3). The final average and SD for each animal was then calculated as the average of the analyzed cycles. The same process was repeated for the program analysis. If the program found a cycle to be "not-representative" via the exclusion algorithm, it was not included in the analysis. Table 3 shows the results between the manual and program analysis as well as the percent difference between the two. Percent differences larger than 10% are bolded. Figure 5 shows the comparison of CFVR per animal between the manual analysis and program analysis. Average HR was equal to or less than 1% different between the two methods. There was a significant difference in both PV and VTI for some animals. Specifically, 3/6 baseline CFP files and 2/6 hyperemic CFP files showed a significant difference in average PV between manual analysis and the program analysis. However, the percent difference between the two methods for these animals Rater was no greater than 17.43%. 4/6 baseline flow patterns and 2/6 hyperemic flow patterns showed a significant difference in average VTI with the largest percent difference of 55.07% found in one animal. The greatest percent difference in the CFVR calculation was 17.02% in one animal while 4/6 animals showed a percent difference of less than 5% between the two methods.

1%-Baseline
Peak Velocity (mm/s) Analysis Duration. Though there were only 6 animals included in this validation study, the manual analysis portion was time and labor intensive mainly because 13 total measurements were made for each analyzed PW Doppler cycle (only three clinically-relevant measures are reported here). On average, it took the raters approximately 1500 total minutes to complete the manual analysis of both viewings. For the program analysis, the raters spent approximately 50 total minutes completing both viewings including time taken to export the raw video files. That is a 30-fold decrease in time spent when using the MATLAB program to measure the specified CFP parameters over manual intervention.

Parameter Values Between Normal and Disease.
We show an appreciable agreement between our program and manual analysis when identifying clinically useful parameters of coronary flow. Additionally, we aimed to determine whether our newly identified parameters varied between normal and disease. Table 4 lists the average and SEM for each parameter as well as the p-value when comparing normal and diabetic mice. There were a number of significant differences especially under hyperemic conditions. All velocity values and HR were significantly reduced in diabetic mice under stress. Decay slopes 1 and 2 were significantly decreased while a significant increase was noted in diastolic decay time 1 for diabetic mice under hyperemic conditions. The only parameter significantly different under normal conditions was HR.

Discussion
The CM is a unique and important microvascular bed that supports the viability of the heart through blood flow and nutrient exchange while constantly experiencing transmural contraction and relaxation forces. In the progression of certain diseases, the CM has been shown to become structurally and/or functionally impaired, often prior to occlusive macrovascular atherosclerosis, which can lead to decreased perfusion, ischemic events, and myocardial infarction. TTDE is a non-invasive clinical tool that has been successfully utilized to assess these impairments of blood flow regulation in the CM. Analysis of coronary TTDE examinations however can be time consuming and prone to operator measurement variability. Automating this analysis could reduce operator variability, allow for the incorporation of more cardiac cycles with additional parameters measured, and reduce labor time. There are, to date, a limited number of studies concerning the automation of TTDE CFP analysis and none, to our knowledge, have been published utilizing an echocardiographic software intended for animal models. We aimed to develop a MATLAB program that could analyze PW Doppler CFP files exported directly from our Visual Sonics Doppler Echocardiography machine system to (1) automatically extract VTI, HR, and newly defined time intervals, velocities, and average slopes from murine CFPs and (2) automatically calculate CFVR. Finally, we hypothesized that this program would significantly reduce analysis time, that measurements obtained from the program would be in good agreement with manual analysis, and that the program would reduce inter and intra-operator variability.
Our results show that our MATLAB program was able to effectively reduce intra-and inter-operator variability when selecting PV, HR, and VTI parameters for analysis. For example, average PV difference between manual viewings was as high as 22 mm/s and difference between raters as high as 13 mm/s. Our MATLAB program, reduced this variability in mean difference to as low as 5 mm/s between manual viewings and 9 mm/s between raters respectively. Reducing variability in TTDE analysis individually and between operators is of critical importance to accurately and precisely make conclusions about blood flow and the CM.
We then investigated how well our MATLAB program agreed with manual analysis in parameter selection by comparing the averages from each animal. While the percent difference was minimal in mean HR for all animals, there were some mean significant differences in PV and VTI. Generally, the program measurements for PV were closer to the manual measurements under hyperemic conditions, with only 2/6 animals being significantly different. Contrarily, 3/6 animals were significantly different for PV when evaluated under baseline conditions. Under 3% isoflurane administration, the coronary arteries dilate resulting in larger blood flow velocity profiles for each cardiac cycle. Thus, the PW Doppler flow patterns often become fuller and more distinct, which is likely why the program better selected average PVs per animal under hyperemic conditions. However, the largest percent difference found for PV in this validation under both 1% and 3% isoflurane was 17.43% (Animal 4, 1% isoflurane), which corresponded to an average difference of approximately 39 mm/s for that animal. The Vevo 2100 software, given the window size that manual PV measurements were made, is accurate to approximately 4-5 mm/s per pixel for the velocity axis. The 39 mm/s difference measured between the program and manual analysis for that animal therefore was only a difference of ~9 pixels. As shown above in the Bland-Altman limits of agreement, the intra-observer variability for PV was larger than 39 mm/s for the raters. Despite the program significantly under-estimating PV in some animals, the deviation was small on a pixel scale and the largest deviation was still found to be less than the intra-observer error.
It is important to note the degree of variability that is present in TTDE recordings, especially with murine coronary flow. Often, the expected biphasic outline observed in coronary TTDE can be partially shadowed making parameter extraction difficult and subjective, even with manual interpretation (Fig. 6). In situations where CFP cycles were partially shadowed, we found that the raters would often extrapolate the CFP envelope to create the expected biphasic outline, basing their outline on flow patterns before and after the cycle in question or from previous experience. This invariably led to the measurement differences observed between our program and the manual analysis.
VTI was also found to be significantly different in some animals. Similar to PV, VTI was more accurately measured by the program under hyperemic conditions because the flow patterns were fuller and more distinct. Under baseline conditions, 5/6 animals VTI was greater than 10% different between the two analysis methods. This was also a result of manual extrapolation leading to an over-estimation of VTI. The largest percent difference was found to be 55.07% which corresponds to an average difference of ~34 mm (Animal 6, 3% isoflurane). While recordings with clearer and more distinct biphasic CFPs would have certainly validated more closely in this animal to animal comparison, CFPs that require extensive extrapolation may not be suitable for the program. To our knowledge, there is no criteria for determining when to exclude a CFP from analysis when the full biphasic pattern is partially shadowed other than at the rater's discretion. Variability would be reduced if discrete guidelines were implemented as to when incomplete CFPs should and should not be analyzed.
We also investigated how well the program measured CFVR. Despite some of the PV values being significantly different, 4/6 animals had differences less than 5% between the two methods. The program significantly underestimated PV for Animal 6 under both baseline and hyperemic conditions, yet the CFVR was nearly equal to the manual assessment. This highlights the importance of consistency when assessing TTDE CFPs. While the program did at times underestimate PV, the automated nature of a program like this will inherently produce consistent measurements across all flow patterns and subjects compared to the subjectivity of individual raters and their potential unconscious bias.
Finally, we introduced the aspect of dissecting and defining CFP cycles further into four phases based on the clear changes in slope during systole and diastole. There appears to be a lack of consistency in CFP nomenclature   14 . Using our program, we found significant differences in time, velocity, and slope parameters between normal and diabetic mice especially under hyperemic conditions. And similar to the aforementioned studies, much of these differences were observed in the diastolic portion of the cardiac cycle. These differences detected by our program may highlight some aspects of early coronary microvascular remodeling that may otherwise go unnoticed in the traditional clinical setting. Future studies will investigate these parameters further to determine if they have any diagnostic or predictive value in determining the onset and progression of microvascular disease in diabetes.
Limitations. We showed that our automated program was able to effectively reduce observer bias, reduce analysis time, and measure PV, HR, VTI, and CFVR with reasonable accuracy compared to manual analysis, especially in clear biphasic TTDE CFPs. There were still significant differences for some animals especially at baseline conditions and when measuring VTI. Future studies should aim to improve upon the algorithm for overlaying the CFP envelopes particularly when certain cycles are shadowed. Potentially incorporating a machine learning algorithm in conjunction with image processing techniques could lead to better envelope overlay and parameter selection.
Our study focused on a single Doppler ultrasound machine used in animal models. Creating a more flexible software that can conform to other ultrasound machine templates would be the next step in creating a comprehensive program. Additionally, there are subtle differences between human and murine CFPs, notably in the systolic portion of the cycle. The algorithm for selecting parameters in this region would likely need to be altered if the program were to be used on human subjects.
Lastly, our program is not entirely automated as the user is still required to manually export the raw video files, input scaling factors, and potentially alter the envelope threshold in the program. Even so, it is much quicker than manual analysis and similar to other non-coronary flow TTDE programs that have reduced time of analysis by 7.5-25 fold 14,20,22 . However, creating a fully automated program would reduce the analysis time even further.

Conclusions
We developed a MATLAB program that can analyze murine PW Doppler CFP video files from the Vevo2100 Doppler ultrasound machine faster than manual analysis while reducing user bias. Building upon the limited number of previously published studies, our program supports raw exported AVI video files which can incorporate a larger number of flow pattern cycles for analysis than single still frame images. We also defined additional parameters of CFPs that may have useful diagnostic implications in the future.